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ABSTRACT 


NUMERICAL SIMULATION 

OF THE NONLINEAR RESPONSE OF COMPOSITE PLATES 
UNDER COMBINED THERMAL AND ACOUSTIC LOADING 

Jayashree Moorthy 1 and Chuh Mei 2 

A time-domain study of the random response of a laminated plate subjected to 
combined acoustic and thermal loads is carried out. The features of this problem also 
include given uniform static inplane forces. The formulation takes into consideration a 
possible initial imperfection in the flatness of the plate. High decibel sound pressure 
levels along with high thermal gradients across thickness drive the plate response into 
nonlinear regimes. This calls for the analysis to use von Karman large deflection strain- 
displacement relationships. A finite element model that combines the von Karman strains 
with the first-order shear deformation plate theory is developed. The development of 
the analytical model can accommodate an anisotropic composite laminate built up of 
uniformly thick layers of orthotropic, linearly elastic laminae. The global system of finite 
element equations is then reduced to a modal system of equations. Numerical simulation 
using a single-step algorithm in the time-domain is then carried out to solve for the modal 
coordinates. Nonlinear algebraic equations within each time-step are solved by the 
Newton-Raphson method. The random gaussian filtered white noise load is generated 
using Monte Carlo simulation. The acoustic pressure distribution over the plate is capable 
of accounting for a grazing incidence wavefront. Numerical results are presented to 
study a variety of cases. 


1 Graduate Research Assistant. 

2 Professor, Department of Aerospace Engineering, Old Dominion University, Norfolk, Virginia 23529-0247. 
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Chapter 1 


INTRODUCTION 


Acoustic fatigue and reliability have been serious considerations in the design of 
aircraft structural components for many years now. Reports of aircraft skin damage 
near jet engine exhausts, prompted the industry to go into an intensive study of 
this problem. Experimental activities were carried out on scaled models as well 
as on prototype assemblies. Analytical studies were performed using simple plate 
models under random pressure fluctuations. Given the limitations of computational 
abilities, these were simplistic linear models of panel response. The predictions of 
panel behavior from these models quite often overestimated response levels. The 
results from analytical work showed very poor agreement with test data. The design 
processes were, therefore, based on empirical relations derived using results from 
tests to modify the simple analytical models. These design guidelines, together with 
improved jet engine configurations, seemed to have alleviated the need for further 
analytical development for some time. 

As modern day aircraft pushed their performance envelope further, the oper- 
ating environment made greater demands on their structures. Apart from large 
pressure fluctuations in their vicinity, aerodynamic heating of exterior skin panels 
implied a high temperature environment. Failure and reliability had to be under- 
stood under combined high thermal and acoustic loading conditions. Moreover, the 
advent of advanced composite laminated structures brought on a new crop of 
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design problems. Their ability to sustain large displacements without yielding, 
moved panel response and behavior into the nonlinear regime. 

All these changes reinvigorated and stepped up new research in the area of 
sonic fatigue. Although, this time around, there was a growing need to develop rea- 
sonably accurate analytical models. With the use of composites, it is impractical as 
well as inefficient to experimentally evaluate every conceivable design configuration. 
Analytical models, as close to being realistic as possible, could be used to numeri- 
cally experiment with a multitude of configurations and predict a potentially good 
design. Towards this goal, a great deal of recent research has tried to refine models 
used for prediction of structural response of panels under random acoustic load- 
ing There are fewer instances of work involving both acoustic and thermal loading. 
The explosive growth in the computational power of present day computers along 
with indepth knowledge in the application of finite element techniques in structural 
mechanics has greatly facilitated this development. 

1.1 Literature Review 

A discussion of the literature on acoustic fatigue must be preceded by a review of 
the research into causes and types of loading. One of the earliest works done in this 
area is a thoeretical study by Lighthill [1], establishing the relationship between 
jet exhaust velocities and intensity of sound radiation. Later on, Lansing et al. 
[2] studied the influence of impinging jet plumes on structures in increasing the 
pressure loads. Another related case occurs when ground reflections of pressure 
waves affect the pressure distribution near a structure. This has been studied by 
Scholton [3]. The second factor giving rise to pressure fluctuations is the presence 
of turbulent boundary layers. The earliest work giving pressure measurements in 
subsonic flow was by Bull [4]. The nature of the pressure spectrum in such a case was 
given by Coe and Chyu [5]. In the same study, the two authors also demonstrated 
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the interaction between boundary layers and high velocity jet exhaust. Pressure 
fluctuations are created due to instabilities such as cavities and shock waves, that 
arise from such mixing. These studies primarily established the sources of acoustic 
excitation in aircraft. Further research has grown into a whole new discipline in 
this area (computational aeroacoustics) and will not be dealt with here. 

1.1.1 Experimental Work 

Investigations into the response of simple structures to loading as discussed above, 
started predominantly as experimental work [6-9]. These, along with many more 
tests were compiled into ‘nomographs’ and design guides [10-12]. Test results on 
honeycomb sandwich panels using fiber reinforced and titanium bonded plates were 
produced by Holehouse [13], and Jacobson and van der Hyde [14]. Composite panels 
were tested by Wolfe and Jacobson [15] , and White [16]. Test data were compared 
with linear predictions in [15], and by White and Mousley [17]. Further experimen- 
tal work [18-20,22], when compared with theoretical results, served to demonstrate 
the shortcomings of linear theories in accurate predictions. This initiated the ne- 
cessity to account for nonlinearities in modelling response behaviors, especially for 
composite structures. Later work [23-29] studied the large deformation aspects more 
thoroughly. A comprehensive review of existing test data to show the presence of 
large deformation was given by Mei [29]. Of these, Bennouna and White [25] com- 
pared nonlinear analytical and experimental results for beams. The paper by van 
der Hyde and Kolb [23] contains test data for a stiffened panel. Mei and Wentz 
[24] compared results from the equivalent linearization method with tests on rect- 
angular panels. WRite and Teh [27] also studied the nonlinear behavior of a panel 
under in-plane compressive loads (less than buckling) and acoustic excitation using 
test results and Rayleigh-Ritz solutions. Robinson et al. [28] also tried to match 
simulation strain spectra with test results for flat and stiffened carbon panels under 


3 



high intensity acoustic loading. This study observed an ‘overprediction in the effect 
of nonlinearity’ in the nonlinear simulation model. The boundary conditions used 
allowed for elastic supports. The pressure time-history used was obtained directly 
from the test readings. 

1.1.2 Analytical Work 
Linear Response: 

The earliest theoretical papers, apart from those with experimental work men- 
tioned above, were all based on linear theories [30-39]. A normal-mode method to 
determine multi-mode response under a space and time correlated pressure distri- 
bution was given by Powell [31]. This was simplified by Clarkson [34]. The studies 
by Samuels and Eringen [32], Bogdanoff and Goldberg [36], Crandall and Yildiz [37] 
and Mercer [38] all dealt with beam structures. The first paper amongst these [32], 
used a generalized Fourier series in combination with a Timoshenko beam model. 
The latter two [36,37] showed that infinite values for mean-square responses under 
infinite mean- squared input existed only for particular models of excitation. Miles 
[30], Eringen [33] and Elishakoff [39] presented results for plates. The Galerkin 
approach was used by Elishakoff [39]. 

Nonlinear Response: 

Among the various methods available for solving nonlinear random vibration 
problems, the approach based on the Fokker-Planck-Kolmogorov (FPK) equation 
was an early entrant. The mathematical background for the application of the FPK 
equation has been given by Caughey [40]. This method generates a governing dif- 
ferential equation for the probability density function of a Markov process. The 
solution to the FPK equation as applied to the Duffing oscillator has been studied 
by Crandall [41], Caughey [42] and Lyon [43,44], among others. Caughey and Ma 
[45] have given exact solutions to a few nonlinear single-degree-of-freedom systems. 
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In the case of continuous nonlinear systems, approximate solutions have been ob- 
tained [46-50]. Nash [46] has compared approximate solutions, using four different 
methods, to a single- degree-of-freedom reduction of the large amplitude vibration of 
plates and shells. He has numerically integrated the reduced equation; represented 
the equation of motion as well as the compatibility equation by a finite-difference 
model and subsequently integrated them; used the FPK exact solution; and finally, 
also used the equivalent linearization approximation. Ahmadi et al. [47] compare 
the FPK approach and the equivalent linearization technique against a perturbation 
method for the nonlinear response of plates under random loads. 

The perturbation technique was extended to stochastic vibration problems by 
Crandall [51]. In addition, multi-degree-of- freedom systems [52,53] and continuous 
systems [46,47,54] have been studied using this method. Another technique for solu- 
tion is the stochastic averaging method, used when there is narrow-band excitation 
of a lightly-damped system [55-58]. A normal-mode approach directly applied to 
stresses has been used by Maymon [59] . 

The one technique that proved to be very effective and convenient in appli- 
cations for nonlinear structured systems is the equivalent linearization technique 
[60-91]. The pioneers in adapting this method to stochastic systems are Booton 
[60] and Caughey [61]. It was then refined and made less cumbersome, for gaussian 
processes in particular [62,63]. Initial application of the equivalent linearization 
method was for single-degree-of-freedom systems [64, 65]; the Duffing oscillator un- 
der broad-band loading [66,67]; and under narrow-band loading [68]. A comparison 
with perturbation solutions for systems with small linearities was given by Crandall 
[64]. Closed form solutions were obtained by Lyon [66], and Iwan and Yang [67]. A 
significant advantage of this method is its easy extension to multi- degree-of-freedom 
systems [69-73]. Spanos [73] and Roberts [74] showed good agreement with results 
from direct simulation solutions. 
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The equivalent linearization technique has been used to study the nonlinear 
random response of beams [75-79] and isotropic plates [46,47,80-84]. These include 
single-mode as well as multiple-mode solutions. Langley [76] presents a formulation 
that directly solves for the response transfer function by using a Cholesky factoriza- 
tion of the equivalent linear equation. Miles [84] modifies the equivalent linearization 
procedure such that the response need no longer be assumed gaussian. He compares 
single-mode results of the fatigue- life for a plate with those obtained from numerical 
simulation and from the probabilistic estimate for damage rate. The analytical so- 
lutions for the response of composite laminated plates under acoustic loading have 
primarily been presented by Mei and his associates [85-91]. All these papers give 
single- mode solutions. The paper by Mei and Wentz [87] studies anti-symmetric 
angle- ply laminates and the rest pertain to symmetric laminates. Of these, Ref. [88] 
deals with nonlinear damping, Ref. [89] studies the effect of transverse shear, Ref. [90] 
uses elastically restrained edges, and Ref.[91] uses mixed boundary conditions. 

Some excellent review papers bring up to date the developments in random 
vibrations until the mid-eighties [92-95]. 

Finite Element Methods: 

Probably the earliest effort to model a random vibration problem using finite 
elements was by Jacobs and Lagerquist [96]. This was followed by a few more 
which were still restricted to linear theory [97-100]. Olson [98] modelled the excita- 
tion cross spectral density over each element with a polynomial function and then 
proceeded with the spatial integration. Dey [99] presents solutions using a matrix 
inversion method and the normal mode method. The global matrices are formed 
using the finite element package ASK A, and subsequently used as inputs in the 
formulation that follows. Apart from these, initial applications of finite elements 
relied on large programs like NASTRAN [18,22,35]. The first instance of a finite 
element model applied to a nonlinear random response problem seems to be in the 
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paper by Hwang and Pi [101]. A conforming triangular element has been used 
to develop all the global matrices for a plate system. The inplane displacement 
equations axe eliminated by static condensation and the equivalent linearization 
technique is then applied. This was followed by work on isotropic beams and plates 
by Mei and his associates [102-105]. All these studies were based on the equivalent 
linearization method. Random response of beams was studied by Mei and Chiang 
[102] and Chiang and Mei [103]. This was followed up by a study of thermally 
buckled beams by Locke and Mei [104]. The coupled thermo- acoustic problem was 
treated in two steps. Thermal buckling analysis solutions were used for the random 
vibration analysis as initial deflections. The study on beams was then extended to 
plates by Locke [105]. The formulations were similar in the two cases. A uniform 
white noise gaussian pressure distribution of increasing decibel levels was used. A 
24-degree-of-freedom rectangular plate element was utilized. Solutions were com- 
pared to their earlier classical single-mode solutions. A convergence study on the 
number of modes active in the solution used up to four symmetric modes. The for- 
mulation used a static reduction or static condensation concept to eliminate inplane 
equations. The mode-shapes for the thermally buckled plate were used to carry out 
a modal reduction of the system equations. 

The next level of advancement was the inclusion of composite laminated plates 
in random vibration analysis [106-108]. Robinson [106] studied the nonlinear ran- 
dom vibration of a particular laminated plate using a time-domain approach. A 
symmetric quarter-plate model of 48-degree-of-freedom rectangular elements was 
subjected to a spatially uniform high-decibel pressure distribution. The nonlinear 
equations were carried through every time-step and a Newton- Raphson iteration 
scheme was adopted to obtain converging solutions within each time-step. This 
simulation, of the nonlinear system of equations, has been implemented using a 
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modified single-step integration algorithm. Chen [107,108] also worked on com- 
posite plates using the equivalent linearization technique. His formulation includes 
consideration of nonwhite gaussian processes and processes with non-zero means. 
He has also presented results for different types of loading assumptions, including 
spatially random gaussian processes. A normal-modes reduction combined with an 
equivalent linearization approach has been verified against solutions from Monte- 
Carlo simulations. 

Numerical Simulations: 

Simulation techniques seemed to have been used popularly as a tool to ver- 
ify analytical results obtained from other methods [46,73,74,84,107,109]. The first 
attempt at formulating a Monte-Carlo technique for a geometrically nonlinear sys- 
tem was likely by Dowell [110]. Most efforts at applying the Monte- Carlo simula- 
tion technique for a random vibration problem concentrate on developing suitable 
methods of load generation. Some of the research in load modelling was not in 
combination with numerical simulation, although they did lay the groundwork for 
further study [33,111,112]. The paper by Dyer [111] models a correlated pressure 
field which is convecting and at the same time decaying over a plate. As the scale 
of correlation distances grows much smaller than the plate dimensions, the pressure 
field reduces to a purely uncorrelated random ‘raindrop’ process as mentioned by 
Eringen [33]. The last of these papers [112] offers another model for the pressure 
load by treating it as a correlated stream of impulses arriving randomly. 

A digitally generated gaussian random load is used by Belz [113] to study the 
nonlinear vibration of a beam. The power-residue method has been used to gener- 
ate a sequence of ‘pseudo-random’ numbers. These are then transformed to a set of 
gaussian numbers, which in turn are scaled suitably to give the random load vector. 
A finite- difference form of the equation of motion is solved assuming that the load is 
a concentrated point load. Robinson [106] uses a similar method for load generation 
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and goes on to use the single-step algorithm developed by Zienkiewicz et al. [114] to 
solve a system of finite element equations. Other than these works, a major contri- 
bution to the field of Monte-Carlo simulation has come from Shinozuka, Vaicaitis, 
et al. The simulation of a multivariate, multidimensional random process using a 
combination of cosine functions with random phase was developed by Shinozuka 
[115-118]. Each of the functions has a weighting amplitude and the random phase 
angles axe distributed uniformly between 0 and 2n. This series of cosine functions 
is summed at constant frequency and wave- number intervals. The weighting am- 
plitude is determined from the cross-spectral density prescribed. The simulation 
process can be applied to non-homogeneous processes, too, the ‘non-homogeneity’ 
being defined by a constant spectral-density that may be made ‘transient’ by a de- 
terministic function of time and space. Another method based on the generation of 
discrete frequency functions corresponding to the desired random process was given 
by Wittig and Sinha [119]. A quick FFT of this generated series then results in the 
actual random process. Sengupta [120] formulated a procedure to derive the spec- 
trum of a decaying, non-homogeneous field. This modelling is particularly directed 
at propeller-driven aircraft excitation, where the excitation signal has dominant 
harmonic components at the blade-passage frequency. 

Solutions based on the Shinozuka approach for various types of panels were 
given by Vaicaitis and others [121-126]. The formulation for the plate governing 
equations is based on an analytical Galerkin method. A fifth-order Runge-Kutta 
integration scheme is then utilized. Results generated cover a broad range of sta- 
tistical characteristics including variances, probability densities, peak probability 
densities, crossing- rates, damage estimates, etc. The papers include response data 
for stiffened panels and composite panels. A variety of loading possibilities due to 
turbulent boundary layers, exhaust noise, unsteady aerodynamic pressure, cavity 
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pressure, shock wave impingement, parametric excitations as well as aerodynamic 
heating have been considered. 

1.1.3 Thermal Effects 

High performance military aircraft of today and the future high speed civil trans- 
port planes can have high acoustic and thermal loads occurring simultaneously 
for significant durations. The effects of high temperature loading on random re- 
sponse, therefore, have become a very relevant concern in aircraft structural design. 
Thermal loading conditions in combination with acoustic loading have been ex- 
perimentally studied before [20,109,127-130], Random vibration experiments were 
conducted on initially buckled beams by Seide and Adami [109]. An attempt was 
made at predicting the onset of the ‘snap- through’ phenomenon. The earlier activi- 
ties [127-130] showed a decrease in the life of all types of panels (isotropic, stiffened, 
composite), at a particular stress level, with increasing temperature. All these tests 
reaffirmed the existence of the ‘snap- through’ phenomenon near critical buckling 
loads and the necessity to develop better analytical models. 

More recently, Ng along with his partners [26,131,132] studied plates under 
thermal [131,132] as well as static in-plane [26] loading. The tests in Ref.[131] seem 
to indicate that severe snap-through and chaotic motion occur more with rigidly 
clamped edge conditions. A little leeway in in-plane motion, as is usually the case 
in everyday structural assemblies, may alleviate such extreme stress cycles. Ng 
and White [26] used the Rayleigh-Ritz method to study a composite plate under 
inplane compressive loads in the buckling/postbuckling range. Comparisons were 
made against test data. Thermal effects on random response were also studied by 
Lee [83]. Apart from temperature variation over an isotropic plate surface, he also 
considers thermal gradients across thickness, or, thermal moments. The analysis is 
a single-mode representation of the equations of motion, solved using the equivalent 
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linearization technique. He notes that the influence of temperature variation on the 
rms displacement and stress values steadily drops off as the intensity of acoustic 
loading increases. In addition, he concludes that although thermal moments make 
little difference in the displacement solutions, normal stresses depend on the tem- 
perature gradient. This becomes significant when considering fatigue life aspects 
for a structure. 

A good overview of research discussed in the above sections is available in a 
report by Clarkson [133] and a recent paper by Vaicaitis [134]. 

1.2 Outline of Present Study 

A Monte Carlo simulation of a laminated plate under simultaneous acoustic and 
thermal loading is the main focus of this study. As mentioned before, the combina- 
tion of temperature loads and high-intensity acoustic pressure fluctuations is a very 
legitimate concern in the development of the high speed transport aircraft. Existing 
literature on this problem still leaves room for good numerical models for prediction 
of structured response under such conditions. Of the very few published works that 
relate to the combined problem, only one [83] considers thermal gradients, using an 
analytical single-mode plate model. None of them utilize the advantages of a finite 
element model. The flexibility in modelling an arbitrarily laminated plate under 
complicated boundary conditions is exclusively a feature of finite element methods. 
As has been demonstrated in earlier research, neither test setups nor real life as- 
semblies show purely simply-supported or purely clamped boundary conditions. A 
certain amount of in-plane motion is allowed. This can be attempted in a straight 
forward manner using finite elements. 

The plate theory to be used accounts for transverse shear rotations and strains. 
This helps application of the formulation irrespective of plate thickness. Further- 
more, accuracy is improved when bending-extension coupling exists, particularly 
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when the laminate has only a few layers. In addition, higher frequencies are pre- 
dicted more accurately as plate flexibility is assumed to be greater in a transverse 
shear model. 

The use of a numerical simulation technique provides a more detailed picture of 
the results than, say, the equivalent linearization technique. A variety of statistical 
information including higher-order moments, probability densities, peak probability 
densities, etc. can be extracted from response time-histories. Moreover, the spec- 
tral densities obtained, more accurately show the widening of peaks expected in 
nonlinear responses. With the computational speed and size of today’s computers, 
an attempt to generate a broad database of results considering various parameters 
relevant to random vibration is a viable project. 

Nevertheless, a modal reduction of the large system of displacement equations 
is carried out. This has the advantage of representing the desired plate system 
with a fewer number of time-dependent equations. A multiple mode analysis for a 
nonlinear system will show the interaction between modes, which is not available 
otherwise. For instance, the relative drop in the rms modal displacement value of 
the first mode in a multi-mode solution, as compared to a single-mode solution, as 
occurs in some cases, can only be observed by multiple mode approximations. A 
further advantage of using a simulation technique is the flexibility available in load 
modelling. A spatially correlated pressure distribution, in the case of a travelling- 
wave acoustic field that is not normally incident, is generated and its effects on the 
response studied. 

As stated above, the temperature distribution is assumed to be more general, 
including gradients across the thickness. The surface distribution is restricted to 
being symmetric, in keeping with data from test setups. The formulation also 
accommodates initial imperfections in the plate els well as initial stresses. The 
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numerical integration scheme used in this study follows the method implemented 
by Robinson [106]. 

An outline of the material in the following pages is now presented. Chap- 
ter 2 develops the finite element formulation of the governing equations and then 
transforms them to a modal system of equations, ready for integration in time. In 
Chapter 3, the single-step algorithm in combination with the Newton-Raphson it- 
eration scheme is developed. The details of generating a gaussian random pressure 
time-history, correlated in one spatial dimension, are explained in Chapter 4. The 
methods used to generate the response statistics and spectrum are briefly explained. 
Other considerations that come up during implementation of the finite element pro- 
gram are also discussed here. These include the determination of mode-shapes for 
systems where in-plane degrees-of-freedom may or may not be coupled with bend- 
ing degrees-of- freedom. Finally, Chapter 5 presents various numerical examples and 
offers inferences. Concluding remarks are presented in Chapter 6. 
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Chapter 2 


FINITE ELEMENT FORMULATION 


Governing differential equations modelling the nonlinear behavior of a plate 
under random dynamic and thermal loads are derived in this chapter. The equations 
of motion axe formulated based on the first-order shear deformation plate theory. 

2.1 Assumptions 

The following are the assumptions find restrictions that go into defining the problem. 

• von Karman large deflections 

• small strains and rotations 

• rotatory inertia included 

• orthotropic, linear elastic lamina of uniform thickness 

• plane stress 

• proportional damping 

The application of the first-order theory, including shear deformation, implies 
that inplane displacements axe linear functions of the transverse coordinate, z. Also, 
straight lines normal to the mid-plane before deformation remain straight but not 
normal af ter deformation. Further, transverse normals axe inextensible, i.e., trans- 
verse normal strain is zero. As part of the plane stress assumption, transverse 
normal stress is treated negligible compaxed to the inplane and the transverse shear 
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stresses. As a consequence of the inplane displacements being linear in z, transverse 
shear strains, and therefore stresses, turn out to be constant through the thickness 
of the plate (hence not satisfying the free surface zero stress boundary conditions). 
To partially rectify this discrepancy, a shear correction factor is used. 

2.2 Strain-Displacement Relationships 

The displacement field over any element of a shear deformable plate is represented 
by 

U(r,y,2r,t) u(x,y,<) z^> x (x,y,t) 

V(x,y,z,t) u(x,y,t) -f* z^y(x,y,t) 

W(x,y,z,<) = w(x,y,t) + w a (x,y), (2.2.1) 

where u,v, w are the plate mid-plane displacements along x,y,z respectively and 
ipxi’fpy are the rotations about y and x axes respectively. u> 0 (x,y) is the initial 
imperfection or the deviation from flatness of the plate. It must be noted here that 
w 0 (x,y) has no restrictions to being smaller than w. 

The above five degrees of freedom, u,v,w,ip x ,ip y , are expressed in terms of 
element shape functions, N } (see Appendix A), i.e., 

u (x, y, t) = Nja"{t) 

r(x,y,t) = Nja^t) 
w(x,y,t) = NjaJ(t) 

0x(x,y,t) = Njaj x (t) 

*l>y (*,y,<) = (2.2.2) 

where a“ , a v - , aj , a^ z , aj v are the values of the displacements u, v, w, tp x , ip y at the 
j th node of the element (see Fig.l). 
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The shape functions, Nj, are chosen to be Lagrangian quadratic interpolation poly- 
nomials [135]. The same shape functions are used for all the five degrees of freedom, 
since isoparametric integration is to be used over each element. 

In general the strain-displacement relationships are written as, 

£,y = + u >.*' + = 1,2,3 (2.2.3) 

where a ‘comma’ followed by a subscript denotes differentiation with respect to that 
coordinate. The notation e* ; represents the Lagrangian nonlinear strain tensor. 

Using deflections in the von Karman sense, it is assumed that components of 
strains due to inplane displacements are much smaller than those corresponding to 
transverse deflection, w(x,y,t). Thus, the expressions for strains reduce to, 

£l = £|l = U„ 

C2 = S22 = V„ + 5 (W„f 


£4 

= 2£23 

= v, z 

+ W, y 


£5 

= 2£3i 

= W, x 

+ U, z 


£6 

= 2£i2 

> 

It 

+ V, X + W, x W,y 

(2.2.4) 


Strains with contracted notation subscripts have been simultaneously defined to 
facilitate further derivation which uses index notation extensively. Substituting for 
U y V, W from Eq. (2.2.1) and writing the shear strain terms separately, we have, 



and 
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Or, in other words, 


e< = t? + ZK{ + e? + e?, 

and 

7m = 7m + 7 mi 


i = 1,2,6 
m = 5,4 (2.2.6) 


where e w is the geometrically nonlinear strain term and the rest are all linear func- 
tions of displacements. Note that in writing strain components from ( w + w Q ) there 
is no contribution from derivatives of w 0 alone since the plate was strain-free with 
initial imperfection, w 0 . Next, the displacement functions above are represented in 
terms of their shape functions. As mentioned earlier, all matrices and tensors (rank 
greater than 2) are described using index subscripts henceforth. Using Eqs. (2.2.2) 
and Eqs. (2.2.5) and Eq.(2.2.6), we get, 


= 

h 

1 t 

0 

Nj, y 
Nj,*. 

«) 

II 

a> 


= CnaP, 



l = j 

(2.2.7) 

Similarly, 


0 

f nf 1 1 



K, = 

0 

Nj,y 

1 1 




Nj,v 

Nj ,x J 

l a } ) 



= C.,a* 




(2.2.8) 


where npe denotes the number of nodes per element which is the maximum value 
j can attain. Since the shape functions Nj are common for all degrees of freedom, 
the matrix Cu is the same for and k,. 



In like fashion we can write 


A = 


Nj tX W 0 ,x 

Nj,y Wo,y 

L^i.y U, °>x Nj,xWo,y] 


{•?} 


C O W 

ij a j > 


(2.2.9) 


C = x < 


[ n,,,n„,' 

Nj,yNk,y 

2 N hX Nk,y 




C w 

ijk a j a k ) 


j-k = l,...,npe (2.2.10) 


7^ = 

/ m 


TVj 0 
0 JV 7 - 


]{$} 


JD'l’ 


m 


= 5,4 (2.2.11) 


and 


U» 

7m 




K) 


= B™ } aj ( 2 . 2 . 12 ) 


2.3 Stress Resultants 

Constitutive relations first need to be developed for individual lamina under plane 
stress, as assumed. Additionally, the lamina stiffnesses are restricted to those for 
an orthotropic lamina. 
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2.3.1 Constitutive Relations 

The therraoelastic stress-strain relationship for an orthotropic material is given by 

’Cu C 12 C*i3 0 0 0 ' £\ aiAT 

Cn C 22 C23 0 0 0 £2 - a 2 A T 

C 13 C 23 C33 0 0 0 I e 3 - aaAT 

0 0 0 C44 0 0 | £4 

0 0 0 0 C55 0 £5 

. 0 0 0 0 0 Ces - ' £6 

These are the lamina stresses given in terms of the strains and the independent 

material constants [C] for any ply in its own principal material directions [136]. 

Using the above equation, the stress-strain relations for a lamina under plane 

stress, result in a reduced stiffness matrix , Qij. Rewriting the constitutive relation 

therefore, we obtain 

’Qu Q 12 0 0 0 I f £1 — ajAT ' 

Q \2 Q 22 0 0 0 £2 - Q 2 AT 

0 0 Q 66 0 0 l £e > (2.3.2) 

0 0 0 Q44 0 I £4 

. 0 0 0 0 Q55 .1 £5 i 

where in terms of the engineering constants, 





Q11 = 

Ql 2 = 
Q 22 = 


Ei 

1 - 1^12^21 
^12^2 
1 - ^12^21 
E 2 

1 - 1 ^ 12^21 


V2\E\ 

1 - 1 / 12^21 


and Q44 = G23, Qss = G31, Qee = G12, remain unchanged. 

For aJaminate of many plies, each having its material axes oriented differently 
with respect to the laminate axes, the constitutive relations must be transformed 
from the principal material directions, 1,2 and 6 , to the global coordinates x,y,z of 
the laminate (Fig. 2 ). Upon transformation, the transformed reduced stiffnesses for 
a given layer are 

[<?] = pt 1 i<?] pr 7 ' 
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where 


[T\- t = [fliPTH*]-' 

with [R] being the Reuter’s matrix as defined by Jones [136]. Transforming the ther- 
mal strain tensor from the principal material directions to the laminate coordinates 
as follows 

a x AT ) 

a v AT > = VrX 1, < a^AT > 

a^ATj \ 0 J 

We are now ready to compute the stress resultants with 

-On O 12 0x6 0 On* c-x — o k AT 

0 l2 0 22 0 26 0 0 — 

0i6 O 26 066 _0 0 < 7 X y — a xy AT 

0 0 0 O 44 O 45 7y* 

- 0 0 0 O45 O55 J k 7 xz 

for the k th layer of the laminate. The expression for strain can be substituted from 
Eqs. (2.2.6) in Eq.(2.3.3) and the final form of the constitutive relations for the k th 
layer is given by 





ijj 

fa 

0 


= Qij ( <" + + <; -«?’*)+ (2.3.4) 

*1 = QL ( 7n + 7n ) (2-3.5) 

with i,j, m,n as defined in section 2.2. Implicit is the condition that the plate 
is strain-free under w a (x,y) and c*. The additional stress term cr* relates to the 
specified initial stress distribution, as had been stated in the problem definition. 

Henceforth all the stresses denoted with index subscripts are in laminate coor- 
dinates, x, y, z. 
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2.3.2 Force and Moment Resultants 

Having written down the expressions for stresses in any ply, the various resultant 
forces and moments that are to be used in formulating the governing equations of 
motion for the laminate are to be defined. These are in units of force or moment 
per unit length , as stresses are integrated over the thickness only. The resultant 
forces, Ni and Q m , and the resultant moments, M, are 


Ni 

M t 


and 





i = 




1 , 2,6 ( 2 . 3 . 6 ) 


(2.3.7) 


(2.3.8) 


where 1,2, 6, 5, 4 now correspond to r, j/,xy, xz,yz directions respectively, and AC 2 
is a constant shear correction factor. Integration over the thickness is carried out 
by summing up the integrations over L layers. The distance from the midplane of 
the plate to the bottom of the k th layer is z*-i, with z 0 = —h/ 2. 

Substituting for stresses from Eqs.(2.3.4-5) into Eqs.(2.3.6-8) above, the stress 
resultants in terms of laminate stiffnesses are as follows. 


II 

£ 

:<? + 


+ 


Mi = Bn | 

:«r + 

*7 

+ 



and the shear force resultants 

m ~ ^ ^mn i 


) + B„n, - N? t + N‘ (2.3.9) 

) + D„k, - Mt T + M; (2.3.10) 

m;n = 5,4 (2.3.11) 
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where 


L r*k _ 

U; = X / Qijdz 

J *h-i 


k = 1 
L 


= X! G?i( z * - z *-i) 


i;i = 1,2,6 


*=1 


axe the laminate extensional stiffnesses, 


“ rzh _ 

Bij = ^ ^ I Qijzdz 
fc=i '' Zfc - 1 

= E*,(^) 

k—1 X 7 


are the laminate bending-stretching coupling stiffnesses, 


z 2 dz 


Dij = E / &*■ 

= E <& (^) 

t=i v 7 


are the laminate bending stiffnesses, 

^ [*k 

N, 


t T = E / 0a4 T ( z ) dz 

k = i 

is the thermal force resultant, and 


Mf T = E / O.j 2 T (z) dz 

*=i • /z *-‘ 

is the thermal moment resultant. .4 m „, the laminate shear stiffnesses, are defined 
identical to extensional stiffnesses, A t j. 

2.4 Hamilton’s Principle 

The governing equations for the problem are derived based on Hamilton s Varia- 
tional Principle, modified by a weighting function, W(<), to be used for numerical 
integration in the time-domain. This can be mathematically expressed as 

& J** W(t) - U - V + Djdt = 0 (2.4.1) 
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where 


T = Kinetic energy of plate 
U = Strain energy of plate 

V = Potential due to external loads on the plate 
D = Dissipative energy of plate; and 

t\ and t 2 imply initial and final time of the time domain. Each of these expressions 
will be derived in the following sections. The dissipated energy is due to structural 
damping, which is assumed directly proportional to the inertia of the plate. 


2.4.1 Kinetic Energy 

The kinetic energy of the plate in terms of the three velocity components can be 
written as 

T = [ - p{U 2 + V 2 + W 2 )dV (2.4.2) 

J\ 2 

where p is the density of the plate and U , V and W are the displacements as defined 
in Eqs(2.2.1). A dot over the variables indicates differentiation with respect to time. 
On substituting from Eqs. (2.2.1) for the overall displacements and carrying out the 
variation operation, we obtain 

fh/2 r . . . 

8T = / p[(u + zip x )(8u + z8rp x ) + (v + zip y )(8v + zSxpy) 

J-h/2 J A 

+ {w8w)]dAdz (2.4.3) 


Note that w Q does not contribute since it is constant in time. All the dependent 
variables in the integrand are independent of z. The following inertias can hence 
be defined 


(/i , h , h) 



z 2 ) dz 


S3 
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As p is uniform for the laminate, I 2 — 0. Next, the kinetic energy term of Eq. (2.4.2) 
is integrated by parts with respect to time using the boundary conditions 

5(U) h = S(U) i7 = 0 
5(V) tl = 8(V) t , = 0 
S(W) h = 5(W) t , = 0 

to obtain 

f T dt = — f f [ ii 8u + I\V 8v + I\W 8w + 8tp x + Ji'Py 8xjj y ] dA dt 
Jti Jtt Ja 

Note that with / 3 not assumed zero rotatory inertia is included as stated in the 
assumptions. Upon substituting for displacements and rotations from Eqs.(2.2.2), 
in terms of nodal degrees of freedom and shape functions, the variation of kinetic 
energy can be written as 

8 [ Tdt = - f MijdjSaidt i-,j = 1,2 ,...,N d f (2.4.4) 

Jti Jti 

where N d f includes the total number of global degrees of freedom for the system 
and Mij is the mass matrix of the finite element model for the whole plate as given 
in Appendix B. It is understood that the area integration over the plate is actually 
the summation of element integrations, or assembly of finite element matrices, over 
all the elements in the finite element model. 

The dissipative energy, D, is now represented in a similar way, as 

8 f 3 Ddt = - f 3 Cij a } 8ai dt (2.4.5) 

Jt x Jti 

where C tJ is proportional to Mij. 

2.4.2 External Work 

This term comprises the work done on the plate by the applied dynamic transverse 
load. If p ( t ) is the random dynamic pressure on the plate, then 

8V = J -p(t) <S(w + w a )dA 
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Carrying out the variation and substituting for w from Eq.(2.2.2), the work done 
term in Hamilton’s principle, Eq. (2.4.1) can be described as 

5 f 2 Vdt = - / Pi(t)Safdt i;j = l,2,...,N d/ (2.4.6) 

Jti Jh 

where Pi is the load vector [see Appendix B]. The details of simulating this random 
load time-history are to be described in Chapter 4. 


2.4.3 Strain Energy 

The strain energy in a plate whose material shows linear Hookean behavior can be 
written in terms of its stresses and strains as 

U = f ^<Tie,dV, * — 1,2, ...,6 (2.4.7) 

Jv 2 

Using the constitutive relation for er; (Eq.2.3.3) and carrying out the variation, we 
have 


SU 


— 8 J 2 [ H" Qmn 7n7ni ] dV 

= / [ Qij^j ^ C » + Q mn 7n h m ] dV 

Jv 

= J [ NiS{t r + tf + e°) + Mi 5ni + Q m 5 lm ]dk (2.4.8) 


where i takes on values 1, 2 and 6, and m takes on values 5 and 4. 

Substituting for the stress resultants JV,-, Mi , Q m from Eqs.(2.3.9-ll) and for 
the strains from Eqs.(2.2.6) the expression for the variation of strain energy now 
becomes 


/ A {[-M«r + «? 

+ [^(<7 + 

+ K 2 A mn ( 7 ^ 


+ ej) + B t jKj - N? t + A 7 ]*(er + 

ej + e?) + DijKj - M* t + M:\8k, 

+ + 7“)}dA 


+ 0 


(2.4.9) 
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■u 


/A 
+ Aij 


( | AijefSe? + Ajj£S<T_ + A tj t°8e? + B tj KjSe? + (N* - N* T )5e"' 
^Sef + AaeJStf + Airf8£ + BijKjK + (N’ - N* T )5e? 

+ AijefSe* + Ajj^W + A^Sef + B ijKj 8e° + (N? - Nf T )St* 

+ BitfSKi + Bj^[8Ki_ + Bitf&Ki + D^Sm + (M- - M* T )8 Ki 
+ K 2 {A mnl thi + A mnl M + A mnl +6 7 - + A mn ^S^)) dA (2.4.10) 


where i,j range over 1,2,6 and m,n each take on values 5,4. The strain energy 
components that are functions of e w , are underlined so as to easily distinguish the 
geometrically nonlinear terms among all the above. Naturally, the one expression 
symmetric in e w forms the cubic nonlinear tensor in the a w nodal displacements. 
The rest of the terms that depend on e w are quadratic functions of a w . 

The next step is to expand all the strains in terms of nodal displacements 
and the corresponding shape function matrices as in Eqs.(2.2.7-12). Once again the 
integration over the complete area of the plate is actually carried out by summing up 
the integrands of every element integration over the number of elements in the finite 
element model. The details of every element stiffness matrix and tensor generated 
as coefficients of nodal displacements are described in Appendix B. 


2.5 Equations of Motion 

With the various energy quantities for Hamilton's Principle now ready, they can 
be put together as follows to obtain the system governing differential equations of 
motion. Eq.(2.4.1) becomes 



Mij'dj + Cijdj 


+ K2ijkio.jO.kai — 


+ (Kij + Kfj - K§ T )a } + K\ l]k aja k 
Pi(t) + P’ - P t AT \saidt = 0 (2.5.1) 


where 
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Mij , Cij axe global mass and damping matrices; 

Kij is the global linear stiffness matrix (including w 0 terms); 

K*j is the global stiffness matrix due to N*\ 

Kfj T is the global stiffness matrix due to AT; 

KUjk is the global quadratic nonlinear tensor; 

K2ijki is the global cubic nonlinear tensor; 

Pi(t) is the global random load vector; 

P? is the global constant vector due to N s (including w 0 terms); and 
pA r is the global constant vector due to AT (including w 0 terms). 

All the above matrices and tensors are of order Ndf. All of them, including those 
derived from the terms in Eq.(2.4.10), are detailed in Appendix B. All, including 
the two nonlinear tensors, are constants independent of displacements. The system 
is now ready for the imposition of kinematic boundary conditions. 

The computational cost in satisfying a weighted average of the above multi- 
degree-of-freedom system over small time intervals until final time, <21 will he as- 
tronomical. In order to circumvent this difficulty, the displacement vector, a , , is 
represented in terms of a series expansion that is then suitably truncated. An ap- 
propriate choice for this truncated series approximation, is a linear combination in 
the mode- shapes of free vibration for the plate. 

The global system mass and linear stiffness matrices are used upon imposition 
of boundary conditions. The mode-shapes required are obtained from the solution 
to 

MijUj -I- Kijcij — 0 i]j = 1,2 ,...,Ndf — Nbc (2.5.2) 

where N^ c is the number of kinematic boundary conditions. This is the governing 
equation for undamped free vibration of the same structural system represented in 
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Eq.(2.5.1). The solution, ay, to this system of equations is assumed to have the 
form 

aj(t) = i = 

Substituting this expression in Eq. (2.5.2), results in an eigenvalue system 

[ Kij — u 2 Mij]aj = 0 (2.5.3) 

that gives an array of eigenvalues, w 2 , which correspond to the natural frequencies 
of vibration for this system. However the amplitudes of vibration, ay are indetermi- 
nate. Information regarding the shape of the response is still obtained by assigning 
one of the amplitudes as a fixed value, usually, unity. A scaled nondimensionalized 
expression of this amplitude vector for each eigenvalue is nothing but the eigenvec- 
tor, also the mode-shape of the response at that eigenfrequency. By determining 
the mode-shape at each eigenfrequency, a square matrix consisting of all the mode- 
shapes can be written as follows 

[$] = [{<£l} {<£ 2 } {4>r} ] 

where {<f> r } is the mode-shape corresponding to the r t}i frequency of vibration and 
is directly proportional to a^. The degrees of freedom eliminated while impos- 
ing boundary conditions are now reinserted to make the number of rows in each 
eigenvector Ndf- By doing this the boundary conditions are accounted for when 
transforming the system in Eqs. (2.5.1). 

The matrix [$] is called the modal transformation matrix for the Ndf x Ndf 
system. This modal transformation matrix now can be used to form the series 
expansion for ay as mentioned before, i.e.; 

ay = <£> 1 (0 + ^>2 92(0 + ••• + <^>r9r(t) + ••• 
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where the modal or normal coordinate, q r , gives the amplitude information at the r th 
mode of vibration. These Eire the new unknowns to solve for in order to determine 
ay. 

The truncation, is in using a combination of a few of these modes, say, 3 to 
6 ( much less than Ndf), to approximate all the Ndf components of aj. Physically, 
such an approximation would be consistent with the nature of vibration response 
of a structure. The dynamic response of a system to any kind of excitation is 
largely concentrated in the neighborhood of the first few natural frequencies of the 
system. Even when the random white noise excitation extends to high frequency 
regions, most of the energy in the response is distributed within the lower frequency 
range. Further, if it is known beforehand that there will be no response at particular 
natural frequencies due to the characteristics of the excitation, that corresponding 
modal coordinate may be dropped from the series approximation. In effect, the 
transformation to modal space lets one tailor the displacement assumption to suit 
the application. The exact considerations involved in choosing the appropriate 
modes to represent the system accurately will be discussed in Chapter 4 under 
implementation. 

Thus, let’s say R mode-shapes are used to express 

a j = 4>j r Qr(t), j = 1, . . . , Ndf] r = R < C Ndf (2.5.4) 

Then 

dj = (f> }r q r (t) (2.5.5) 

and 

aj = 4>jr Qr(t) (2.5.6) 
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To obtain the reduced modal equations, the above transformations are substi- 
tuted for the displacement, velocity and acceleration in Eq.(2.5.1), such that 


Mijdj = q r , (2.5.7) 

Cijaj = Cij<t>jr Qr , (2.5.8) 

( K,j + K-j - Kfj T ) a, — ( Kij + K'j — Kfj T ) 4>jr <?r , (2.5.9) 

Klijk ajOk = Klijk<f>jr<i>k,qrq t , and (2.5.10) 

K2{jkiajdkai = K2iju4>jr4k.4>itqrq.qt (2.5.11) 


where r, s,t now range from 1,. In addition, to reduce the matrices in their 

i coordinate (which still goes from 1 ,...,Ndf), the whole equation in Eq. (2.5.1) is 
premultiplied by the transpose of [$]. If <f>i p is used to denote the transpose (index 
i must be used this way, to be consistent with the index in all the system matrices), 
then Eq.(2.5.1), upon substitution of Eqs. (2.5.7-11) and multiplication, becomes 

JV(t) ^ \Mij4>)r4 > ip\qT + [ C,j(f)j r 4>,p j (jr 

+ [(A^ + A ij — 4>j r (f>ip] q r 
-I- [ K2,jki<f>j r <t>ka4>it4>ip] qrq 3 qt 

- «(()*, + P'^p ~ P. AT M dl - 0 (2.5.12) 

It must be clarified that unlike in the case of the familiar matrix notation it is 
not required for 4>i p to be placed before the matrices it ‘premultiplies’. All that is 
essential in the index notation is for the index assignments to be consistent so that 
the right elements of two matrices multiply each other. 
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Therefore, the final form of Eq.(2.5.1) transformed, and thereby reduced, from 
the global system physical degrees of freedom to the new modal degrees of freedom- 
can be now written as 


w{t) | m pT q r + c pr q r + ( k pr + k’ r - k^)q r 


+ kl pr3 q r q 3 + k2 pr3 tqrq»qt - Pp(i) + Pp ~ P 


AT 


| dt = 


0 (2.5.13) 


where p, r, s, t range from 1, . . . , R. 

The matrices m pr , c pr , k pr , . . ., etc. are the modal equivalents of the global 
system matrices as shown in Eq.(2.5.12) above. In actual implementation the modal 
transformations are not done at the global level. There is no possibility of storing an 
Ndf x Ndf x Ndf x Ndf cubic tensor on most machines and neither is it efficient. Instead, 
transformations are carried out at a much smaller element individual degrees-of- 
freedom-level. Details of the implementation process will be discussed in Chapter 
4. 

This is the set of equations that will be integrated from initial to final time. 
The solution vector qr{t) obtained is then used to approximate the displacement 
vector, Oj, as in Eq.(2.5.4). 


w 
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Chapter 3 

NUMERICAL INTEGRATION 


The integration technique used to march forward in time is described here. A 
single step algorithm, as put forth by Zienkiewicz et al. [114], has been adapted to 
generate a system of nonlinear algebraic equations that are solved iteratively within 
each time-step. The salient features of this procedure are detailed in the following 
sections. 

3.1 Single Step Algorithm 

Within a single time-step, say, between time f„ and f n + i, the modal coordinate, q , 
is expressed as a polynomial series of degree p, in r, as follows 


y , a <p) tL 

?( r ) - Xj <f r m + D l 


m=0 


with 


f„+i = t n + At and 0 < r < At 


Or, on expansion, 


q( T ) = 


„ T 2 (p) tP 

q n + <jn T + q n 2J- + • ■ • + a " ~pi 


(3.1.1) 


(3.1.2) 


the dot denoting differentiation with respect to time. This series expression for q(r) 
is exact. 


34 


For a second order system of differential equations, with initial displacement 
and velocity being the natural initial conditions, a cubic degree polynomial is 
the lowest applicable. Any polynomial of higher degree requires initial conditions 
for higher derivatives of q(r) and does not have unconditional stability charac- 
teristics even for a linear system. So, when degree of polynomial in Eq.(3.1.2), 
P — 3, 

•• r2 (3) T 3 

9n+l = <7n 4" Qn T "I" Qn 4" On ~^T 

—2 

. . „ ( 3 ) T 

Qn+l — 9n + QnT + O n v — 

q'n+1 = q n + ct n (3) r (3.1.3) 

When t = At, this becomes 


9n+l = 


9n+l = 


9n+ 1 = 


Qn "h Qn Qn 

Qn + a n 3 -g- 

Qn "1“ Qn Af H" Q, 

5 + Q < 3 > — 

t/n i a n 2 

9n + a n (3) At 
4n + On (3) At 


2! 


( 3 ) 


At 2 


(3.1.4) 


which can be generalized as 

{q}n+l = {q}n + {<*}n 


where the right-hand sides above consist of a known function from the previous time- 
step, {q} n , analogous to a ‘predictor’, and the unknown {a}„, the ‘corrector’. This 
corrector, with the only unknown, a n ^ 3 \ is determined by a weighted satisfaction 
of the modal system of equations, Eqs.(2.5.13), over every time-step. The weighting 
function, W, is defined such that 
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(3.1.5) 


At 

/ W(r)T m dr 
® = e m At m 

f W{r)dT 

o 

where 0 m are free parameters chosen based on stability criteria. These parameters 
can correspond to any of the various numerical integration algorithms such as the 
Houbolt algorithm, the Bossack-Newmark algorithm, the Wilson-0 algorithm, etc. 
In this work the 0 m values were chosen based on the Wilson-0 integration algorithm 
[137]. For a second order linear system, the cubic algorithm (p = 3) needs three free 
parameters, 0i, 02, $ 3 , which for Wilson’s algorithm, take on the following values 

01 = Ow 

02 = 0vv 2 

03 = 9w\ Ow > 1-366 (3-1-6) 

where 6 W corresponds to the Wilson-0 parameter. The convenient fashion in which 
the three free parameters above are defined in terms of the Wilson-0 parameter is 
the reason this algorithm was chosen over the others. This is particularly useful for 
the nonlinear system we have, where the quadratic and cubic terms in q(t) when 
expressed in terms of Eqs.(3.1.3), will result in more than three free parameters 
(see Appendix C). Thus for the nonlinear system of equations, the definitions in 
Eqs.(3.1.6) are extended as follows 

0 m = 0w m , m = l,2,...,3p and 6 0 = 1 (3.1.7) 

In our case 3p is 9. 

The only other function of time in the weighted integral equation (2.5.13) is the 
dynamic load, p(t), which needs to be interpolated appropriately within the time 
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interval, At. It is reasonable to approximate that it varies linearly between f n and 
tn +1 to give 

p„(f) = 9\ p(t„+i) + (1 — 9i)p(t n ) (3.1.8) 

The system of modal equations in Eqs.(2.5.13) is now rewritten with Eqs.(3.1.3) 
substituted for the unknown modal coordinates, q(t). The resulting expression and 
its integration is detailed in Appendix C. As shown in Appendix C, after integration 
has been carried out, the quantity a ( „ 3) is replaced by a new variable defined as 



This has been done for the sake of implementation. For a random dynamic problem, 
q ( „ 3) , analogous to the rate of change of acceleration, is a very large quantity and 
At is usually a small quantity. By solving for d n we are dealing with quantities 
on the order of displacements and are directly obtaining the increments in modal 
displacement, q(f). 

Thus the numerical integration algorithm used herein results in a set of coupled, 
nonlinear algebraic equations in d as given below. 

A2i jkl djd k d, + Alijkdjdk + Aij dj + AOi = *i(d) = 0 (3.1.10) 


where 

A2,jki - arises out of the cubic modal tensor, k2iju 

Alijt - arises out of quadratic functions of dj in k2ijki and the quadratic modal 
tensor, kl(j k 

Aij - arises out of linear functions of dj in k2{j k i and fc l*>Jfc ; the linear modal 
stiffness matrices, k i} , fc/y, fc$ T ; and the mass and damping matrices, M,y and C x j 

AOi - arises out of all the constant terms from all the above tensors and matrices 
that are functions of solutions from the previous time-step, {q}, as well as from the 
modal load vectors at the current time-step, p,(<), and p\. 
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3.2 Iteration Procedure 

The solution procedure for the nonlinear equations obtained above can be chosen to 
suit the degree of nonlinearity of the problem. When a system is weakly nonlinear it 
may be sufficient to solve a linearized set of equations at every time-step, using the 
solution from the previous time-step. In this work, however, nonlinear equations are 
retained, so as not to restrict the scope of this formulation to problems with small 
nonlinearities. Accordingly, the well-known Newton-Raphson iterative scheme has 
been utilized within each time-step. 

If it is required that Eqs.(3.1.10) be satisfied at the (p + 1) iteration, then we 
can write 

*? +1 (d) = 0 

which can be expanded as a Taylor’s series about , from the p t ^ 1 iteration as follows 


*? +1 (d) = tff(d) + Ad r + ••• = 0 ( 321 ) 


By truncating this series to second-order errors, the above is rewritten as 


f)\ I/P 

» - *f (d) 


- ; 

i-e., 


= : 



A 

KT(<F)£ur r = 

(3.2.2) 


where 


i i 

d? +1 =dP + Ad? 

(3.2.3) 

a i 


B ! 

and Kf T is the tangent stiffness matrix. 


= ! 
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It is obtained by substitution of Eqs. (3. 1.10) and differentiation (as in Appendix 
B). i.e., ’ 

= [ A2ijki d p d p d p + Al ijk d p d p + Aijd] J + A0,(q)] Ad p 


A2>irkl "h A%Hrk 


+ A2 iklr ]d p d p 


+ 


A\ir k " 1 “ A\i k r 


d p + At 


|a4 


(3.2.4) 


The corresponding right-hand side vector, — 'J'J’, is computed using d? r and the 
new system of equations in Eqs. (3. 2. 2) is iteratively solved for A d?. The value of 
is now updated using Eqs. (3.2.3). The iterations are carried out until a certain 
convergence level is met. Convergence is said to have been achieved when the error 
between solutions of two successive iterations is within some predetermined bound. 
The error norm chosen here is the maximum norm, i.e., 


jP+i 


max 


d p 


d p+1 


< err 


(3.2.5) 


where err is a sufficiently small positive number. 

The startup of the iterations needs a suitable guess value, which in this case 
can either- be zero (giving the linear solution at every time-step) or the modal 
displacement from the previous time-step, q n . 
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Chapter 4 


IMPLEMENTATION AND SOLUTION PROCEDURE 

The previous two chapters have described how the plate vibration problem 
has been modelled. In order to proceed with solutions to specific problems various 
preparatory computations need to be carried out. These include solving eigensys- 
tems to determine the fundamental frequency, the modal transformation matrix as 
well as the critical buckling temperature all of which axe subsequently required. 
Apart from these, time-histories of random gaussian acoustic pressure distributions 
at different decibel levels need to be generated. Furthermore, post-processing of the 
resultant displacement and strain time-histories requires computation of the power 
spectra, the probability distributions and also the various statistical parameters 
such as mean, variance and other higher order moments. This chapter deals with 
these above-mentioned aspects and certain implementational considerations that 
need to be highlighted. The load generation details will be discussed first. 

4.1 Random Load Generation 

Since we are dealing with a random stationary process , a Monte Carlo simulation 
is used to generate the random pressure time-history. A Monte Carlo simulation 
evaluates the required function at a number of sample points all lying within some 
domain of interest. These values are then used to represent a close approximation 
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of the actual function itself. The error in the estimate of the function is greatly 
dependent on the number of samples used. 

The pressure time-history is a filtered normal distribution whose sound pressure 
level (SPL) is suitably scaled. A uniform random array of numbers is first created. 
This uniform array is then transformed to a normal distribution of numbers with 
zero mean and variance unity. The method used is the KR algorithm as given by 
Kinderman and Ramage [138]. 

4.1.1 Filtered Random Process 

The array of normally distributed numbers, Y , is passed through a recursive linear 
filter. The filtering process is based on the bilinear transformation method [139]. 
This technique maps the frequencies, /, into a new variable, w, using the equation 

w = tan [n(fAt)] (4-1-1) 

where At is the sampling interval or time-step. This results in a frequency response 
function for the filter as follows 

w>i J = (m) (=n^) (41 - 2) 

where w is as defined in Eq.(4.1.1), and a and b are the values of w corresponding 
to the lower cut-off frequency, //, and the upper cut-off frequency, respectively. 
The filtered process can be represented as 

m n 

X„ = y>Y„_ k + Yd,x„- j (4.1.4) 

Jt=o i=i 

where there are m -j- 1 coefficients cjt and n coefficients dj that are constant for the 
filter and define the relationship between the broadband input, Yk and the filtered 
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output X n . These constants are determined using the poles of Eq.(4.1.2) that result 
in a stable H(f). They are 

b 

°° ~ (1 +«)(!+&). 


c\ = 0 


b 

~ (l + a)(l+6) 


(l+q)(l-fr) + (1 — fl)(l + b) 
(1 + o)(l + b) 


, _ ( 1 -«)(!-&) 
2 ~ ( 1 + «)(! + *) 


(4.1.5) 


With these constants and the input array, Y, we can generate the output 
filtered, band-limited process, X, using Eq.(4.1.4). This filter as such does not 
produce sharp cut-offs. The filtering process can be repeated until the desired 
drop-off is achieved. 

The filtered process, X, now needs to be scaled to a gaussian pressure time- 
history. In order to do this the root mean square value of the pressure distribution, 
Prma, needs to be computed. If the level of the acoustic pressure, supplied in 
decibels, dB v is SPL, then 

PL = Pl., 10 Sf>t/ ' 0 ( 4 . 1 . 6 ) 

is the magnitude of the power spectrum levels for the broadband gaussian acoustic 
pressure distribution, in psi 2 /hz. It is a constant through the frequency spectrum. 
The reference pressure level is p re /. The expression for p rms is therefore given by 

pL. = t m!)\ 2 "‘pl,d! (4.1.7) 

Jfl 
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where nf is the number of filters used. The limits // and f u are as defined before. 
When nf is 1, the expression reduces to the familiar mean-square relationship be- 
tween any input and output process. As mentioned earlier, the filtering process 
needs to be repeated in order to produce a more ‘ideal’ band-limited output which 
has a steep drop-off at / u . This repetition is equivalent to applying the magnitude of 
the frequency response function, \H(f)\ 2 , nf number of times. Hence the equation 
for mean-square of the output process is as given in Eq.(4.1.7) above. 

The pressure time-history vector is thus obtained by 

p = p rma X + p (4.1.8) 

where p is the bias or mean, if any. The notation p stands for an array of acoustic 
pressure values. It represents 

T T 

P = {p(<l) P(*2) P{*n)} = {Pi P2 PN } 

which is the expression for pressure used in deriving Eq.(2.4.6). The subscript N 
here denotes the final time-step in the time-marching scheme. 

4.1.2 Normal Random Load 

The random pressure time-history is now ready. The load on the plate must now 
be computed. This load is modelled in two ways. The simpler of the two, due to 
waves along the 2 -axis, perpendicular to the surface of the plate, is treated as being 
due to uniform pressure over the whole surface area at any given time. This model 
will be dealt with first. 

The global load vector corresponding to unit pressure, {F } , is defined in Eq.(2), 
Appendix B. It will be used here as Fi , where subscript i ranges over the global 
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degree of freedom numbers. As shown in the same equation, Eq.(2), the final form 
of the global load vector is given by 

p«(t) = p n (t) Fi , i = 1 , 2 , ... ,Ny, and n = 1 , 2 ,..., IV (4.1.9) 

It must be made clear that N d j is the total number of degrees of freedom in the 
global system and N is the total number of time-steps. This is the load vector, 
Pi(t), used in Eq.(2.5.1). It is obtained simply by scaling Fi by the p rms value from 
the vector p for the corresponding time-step. 

The entire formulation so fax has been based on this load model. At a particular 
time-step, a single p r ms value is used over the whole plate. But if the waves were 
incident at an angle to the plate, say tangential to the surface, then this would no 
longer be an accurate model. 

4.1.3 Space-Time Correlation 

A pressure field at an angle to the plate is treated as a series of plane waves incident 
on the plate at a certain angle, /?, as shown in Fig.3. It is assumed that the waves 
travel from left to right over the plate. Therefore the spatial direction of interest is 
along the x-axis only. The pressure distribution along the y-direction is restricted 
to being uniform. 

The above plane wave acoustic pressure field satisfies the one-dimensional wave 
equation. The d’Alembert solution to the wave equation results in an expression 
for the pressure, p , at any point as follows (outgoing wave term only) 

p = p{t - -) (4-1.10) 

c 

where s is the distance measured along wave direction ; c is the speed of sound in 
air. 
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In terms of the plate coordinate x, this becomes 


, x sinB , 

p = p{t — ~r-) 


(4.1.11) 


This expression essentially correlates pressure values along the x-axis with time. 
The time-lag in travelling a distance x is given by xsinP/c. In other words, the 
magnitude of pressure at the ‘leading edge’ is the same as the magnitude at a 
distance x from the leading edge after a time interval xsinpfc. 

Thus in order to create the global correlated pressure vector over the surface 
of the plate, we first define a global array, g k , as follows 

I* sinp/c 


g k = integer 


At 


k = 1,2 ,...,Nn 


(4.1.12) 


where Ns is the number of nodes in the mesh (Ndf = 5 Nn)- This array contains 
the integer values of the time-lag in multiples of the time-step, At. For example, 

let 

5i = {0 0 1 2 2 3 4 ...} T 

correspond to global nodal x-coordinate values 


x k = Xi, x 2 , z 3> x 4 , x 5 , x 6 , , xi = O.U 

In this instance, node 3 at a distance x 3 from node 1 has a corresponding g 3 value 
of 1. This implies that the time-lag in a wave travelling from node 1 to node 3 is 
approximately equal to At. Similarly, a wave reaches node 4 after a time span of 
2At seconds. 

With this array, g k , defined, the load vector at the n th time-step can now be 
determined. Using the same notation as in section 4.1.2, it can be written as 


P"(t) = p n -g k (t)Fi , 


1,2,..., Ndf, 


k = 1,2,..., N n (4.1.13) 
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with Pn-j* = 0 for all n — </* < 0 (negative subscripts). 

Note that the subscript i represents the global degree of freedom number 
whereas k stands for the global node number. It is understood that gk values 
remain unchanged for all the five degrees of freedom at one particular node. Again 
the example from above is used to explain Eq.(4.1.13). For the degrees of freedom 
corresponding to the first four nodes, whose x-coordinates are x \ , x ^ , X3 , X4 , and gk 
values are 0,0, 1,2 respectively, the final form of P, at the n th time-step is 


Pi ) 


0 'j 

P 2 


0 

p 3 


p n p 3 

P 4 


0 

P 5 


0 

0 

; 


0 



PnPs 

0 

0 

. 

\ — i 




0 

0 

Pn-lPl 3 



0 



0 

0 

0 

Pl8 


Pn-2pL8 

Pi 9 


0 

‘ P 20 i 

1 

0 


The elements of P, and P, axe grouped as 5 degrees of freedom per node, with the 
entries corresponding to the a w degree of freedom being the only non-zero entries, 
i.e., F w ^ 0, and F u = F v = F^ = F*’ = 0. 

Thus, Eq.(4.1.13) gives the expression for the correlated random load vector. 
Had the angle of incidence, ( 3 , been zero, implying normal incidence, all gk values 
would have been zero, resulting in Eq.(4.1.9) for Pi(t). 
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4.2 Implementation 

As mentioned earlier, the preparatory computations that need to be carried out in- 
clude the determination of the modal transformation matrix and the corresponding 
eigenvalues. The appropriate ‘modes’ or eigenvectors are extracted from this matrix 
to represent the relevant degrees of freedom in the response of the plate system. 

It is extremely important that care be exercised in choosing the appropriate 
modes to represent the structural response adequately. There are particularly sig- 
nificant factors which influence the choice of mode-shapes during implementation. 
Important considerations are the damping, spacing between modes and the fre- 
quency range of interest. Also to be included are the properties of the plate, the 
type of dynamic loading and the nature of the stiffness matrix, The last will 

be discussed in detail in the following sub-sections. 

The geometry of the plate and the boundary conditions on its edges affect the 
frequency content in the vibration response of the plate. There are two planes of 
symmetry of concern, namely, the x — z and y — z planes. If the plate were to be a 
simple rectangular symmetric laminate with symmetric boundary conditions, then 
it is symmetric about both these planes. Its response to a normal random load 
would be restricted to modes that are symmetric in xandy-coordinates. In the case 
of a load that is correlated along the r-coordinate as described in section 4.1.3, the 
only plane of symmetry is the x - z plane even for a fully symmetric plate. This 
still permits the use of a half-plate model, and correspondingly, mode-shapes that 
are symmetric in the y-coordinate. For the case of a plate that has no plane of 
symmetry as in the case of mixed boundary conditions or non- rectangular plates, 
no assumptions in choosing mode-shapes can be made. All the predominant lower 
modes may need to be included. 
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4.2.1 Uncoupled Linear Stiffness Matrix 

In Eq.(30) of Appendix B, the element linear stiffness matrix, [ K e ], is a fully loaded 
matrix only if bending-extension coupling, arising due to matrix in the laminate 
stiffnesses, and initial imperfection, w 0 (x, y), are non-zero. In the case both are zero, 


m 


k mm 0 0 

0 ks ww ks w * 

0 ks k** + ks ** 


where the membrane displacement sub-matrix, [fc mm ], is completely independent 
of the other sub-matrices that correspond to bending and curvature. In this case 
when solving Eq. (2.5.2) 


Mij’dj + Kijdj = 0, i\j = 1,2, - N bc 


at the global level, to obtain the transformation matrix, [$], the eigenvectors cor- 
responding to inplane displacements, {a m }, will be independent of those corre- 
sponding to bending and curvature, {a"'}, {a^}. Thus it is important that {a m } be 
represented by an adequate number of its eigenvectors in the truncated series for 
Eq.(2.5.4). 


In effect, for the case of an uncoupled stiffness matrix, Eq.(2.5.2) splits to two 
smaller systems as follows 


[M mm ]{ a m } + [A' mm ]{a m ) = 0 (4.2.1) 


and 


M ww 

0 


0 

M++ 


{ 5 } 


r Ks ww Ks w + ] f a w \ 

+ [ Ks' f,w K++ + Ks** J \ a* J 


= 0 (4.2.2) 


The global sub-matrices used in this equation axe used as direct assembled exten- 
sions of the element sub-matrices defined in Appendix B. Solutions of these two 
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eigenvalue systems result in expressions for the nodal degrees of freedom as given 
below 

{a m } = ( 4 - 2 - 3 ) 

= [*^>]{ 9 ( *^} (4-2.4) 

where [$ m ] is the modal matrix for the membrane system, Eq.(4.2.1), and {q m } are 
the corresponding membrane modal coordinates. Similarly, [$(“’^4] is the modal 
matrix for the bending system (o'", o^) and {q( w ^} its modal coordinates. The 
equations above imply that the following series expressions can each be truncated 
to represent the membrane and bending degrees of freedom, respectively. 


{a"*} = {<f> m )iq? + {<T}2<Z 2 m + {*~}nT + ... 

and 


j£\ = + {*<-*>}* df*’ + + ■■■ 


These two equations when put together in matrix form result in 


f {« m } ) 




0 


{ i. 




(4.2.5) 


which is the same as 

{a} = [{4>}\ {<£} 2 {<^}3 {^}4 ... ] 4 


91 

92 
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The length of each mode-shape vector, {</>}, is Ndf ~ N\, c . They are brought back 
up to length Ndf by inserting the kinematic boundary conditions at the appropri- 
ate degrees of freedom. Then the transformation is similar in form to Eq. (2.5.4) 
repeated here for reference 

a i — j = 1,2 ,...,Ndf, r = 1,2 , R < Ndf 

The actual number of modes used, R , is decided based on convergence studies 
for both the menbrane and bending modes, done separately. The minimum number 
of membrane modes required for each problem is determined based on the con- 
vergence of solutions with increasing number of modes. Once that is fixed, those 
bending modes over which a significant response level is observed are retained and 
the rest discarded. This is again done by checking the convergence of solutions. 

For plates with stiffness matrices as discussed above (symmetric laminates), 
there are associated buckling problems. A buckling eigenvalue problem is solved to 
determine the critical buckling temperatures using the equation 

Kij - - 0 (4.2.6) 

where a nominal value of AT is used to compute Kfj T . The critical buckling 
temperature is determined by scaling the AT value used, with the lowest eigenvalue. 
These buckling temperatures are used to evaluate the response of a plate under 
thermal and acoustic loads around these critical values. 


4.2.2 Coupled Linear Stiffness Matrix 


In the case of a laminate with non-zero bending-extension coupling ( Bij ^ 0) 
and/or with a prescribed initial imperfection, w 0 (x,y'), the element linear stiffness 


matrix, [K e ], has the form 




k om 


k *m 


k m0 

k°° + ks ww 

+ ks^ 


k m <P 

k + ks w ' 1 ’ 
■k** + ks++ 
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as given in Appendix B. The sub-matrices in the membrane positions, [fc mo ] and 
[k om ] arise due to w 0 (x,y), and and [fc^ m ] arise due to B i} . Either way, the 

membrane components in [K e ] are no longer independent of the bending compo- 
nents. In this case, therefore, the complete coupled system given in Eq.(2.5.2) is 
solved to give 


1 


,V> 


) = [{ih {4>h {<^>3 


{*} 4 •••] 


9i 

Q2 

Q3 


(4.2.7) 


which is in the same form as Eq.(2.5.4). The required number of modes are used 
again based on modal convergence studies of solutions to the random vibration 

problem. 


4.2.3 Assembly of Modal Matrices 

Once the required modal transformation matrix is ready, all the matrices and tensors 
in the random vibration problem need to be reduced to their modal equivalents. As 
stated at the end of chapter 2, it is physically impossible to assemble K\, jk or 
K2{jkl to the global nodal degrees of freedom level. In fact it is impractical even at 
the element assembly level to store, say, Kl' jk , as a 45 x 45 x 45 tensor. Therefore, 
at the programming stage, the modal matrices and tensors are all prepared at the 
element sub-matrix level. For instance, [fc2“*], a9x9x9x9 tensor for the a w 
degree of freedom is transformed to its modal counterpart as follows 


[k2 ww ] modal = [<t> w '} T [k2 ww ] [<r'H<r e ][<r'] 


and, likewise 


= [ ( f ) "'') T [kl mW }[<i> We ][<f> We ] 


(4.2.8) 


52 


Ml Ell III 


where [<£*"*] and [<£"*'] contain the element modal transformation matrix components 
for the a w and a m degrees of freedom, respectively. All the other element sub- 
matrices defined in Appendix B, are transformed to their modal equivalents in a 
similar fashion. Once all the stiffness, mass and thermal matrices and tensors are 
in modal form, they are then assembled at element level. Since they are no longer 
coefficients of ‘nodal’ degrees of freedom, the ‘assembly’ is simply a straight forward 
summation. The same procedure follows when summing up modal contributions 
from every element to give the ‘global’ modal matrices and tensors. 

The advantage in using modal transformations is evident right here. All the 
way from the sub-matrices level we are dealing with arrays whose maximum size is 
set by the number of modal coordinates used. This smaller system is going to be 
integrated forward in time. The variables of interest, including strains and stresses, 
can be written out at any instant of time desired. 

Assembly of Load Vector 

As described in section 4.1.3 of this chapter, the exception to the above trans- 
formation and assembly procedure, is the formation of the correlated load vector. 
In order to account for the space-time correlation, the global array assembled from 
all elements, F,, is still in ‘nodal coordinates’, calculated for unit pressure per unit 
area. This vector is then prepared for each time-step using the method described. 
Then, a modal transformation is carried out at every time-step. 

4.2.4 Step-Size Considerations 

One of the important aspects of numerical simulations is the time-step size or sam- 
pling interval, At. It is required that At be small enough to capture the response 
up to the upper frequency cut-off, f c . This cut-off frequency or Nyquist frequency is 
the “highest frequency that can be reproduced from data sampled at equal intervals 
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Af” [140]. Therefore, in order to avoid aliasing (where two disparate frequencies 
become indistinguishable), a sampling interval of 

At < ^ (4-2.1) 

A/c 

must be used. 

A good estimate of the f c value to be used can be obtained from the natural 
frequencies of the plate. It is generally known that the major response of a vibrating 
system falls within about its first six modes of vibration (six symmetric modes in 
case of a symmetric system). If the sixth natural frequency is known, then an 
f c value greater than this can be used. The sampling interval can be determined 
accordingly. The load spectrum also is set to span all these frequencies. 

A greater constraint on the step-size for numerical integration is its stability. 
Although there axe no guarantees for a nonlinear system, an adequately small value 
has to be used. This usually results in a sampling interval smaller them that needed 
to accommodate the Nyquist frequency. A point to mention in this context refers 
to some work on stability of nonlinear systems. There are some proofs for stabil- 
ity of nonlinear structural dynamics based on the principle of energy conservation 
[141,142], 

One more note regarding the step-size refers to the membrane plate responses. 
The inplane degrees of freedom for a thin plate are typically active at very high 
frequencies. A At value that captures such high frequencies (at modes that couple 
in case of a nonlinear problem) will have to be extremely small. This results in 
a computationally enormous task. Since the response at this range is negligibly 
small, it is not critical if the vibrations at these frequencies are ignored. All we 
really need, again for a nonlinear problem, is the effect of membrane coupling on 
the lower bending modes. Hence, the step-size required is computed based on the 
significant bending modes only. 
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A change of variable brought about by implementation considerations is the 
replacement of a ^ by d n in the expressions for acceleration, velocity and displace- 
ment at any time-step, see Eqs.(3.1.4). Here, 


d n 


a 


O) 


At 3 

6 


As mentioned earlier, oL 3) is usually very large and At 3 is very small. This change 
of variable avoids the possibility of errors induced by computations involving such 
large number of significant digits in both these quantities. 


4.3 Postprocessing 

The large number of samples that form the time-history data sets for each output 
variable in the solution need to be processed to extract meaningful information re- 
garding plate response characteristics. Each of the parameters used will be discussed 
below. 

The power spectrum estimate, P(f), is based on the description in [139]. 
For a data set, Cj(t), of N samples, the FFT using a ‘Parzen window’ function, is 


N-i 


Ck = J2 c > w ) 


,2nijk /N 


k = 0 - 1 


where 


i=o 


Wj = 1 — 


3 - i(N - !) 

} (AT + 1) 


(4.3.1) 


(4.3.2) 


is the window function. The power spectrum estimate, therefore, is 


p{fk) = ^ (|c*l 2 + |Cn-*| 2 ) , 


k = 1,2,..., (y 


where 


N 


W = N w ) 
>=o 


1) (4.3.3) 


(4.3.4) 
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Frequencies /* are defined only for zero and positive frequencies as 


fk = 


N&t ’ 


N 

k — 0, 1, ■ ■ • j g 


(4.3.5) 


The power spectrum is finally obtained by averaging together K estimates at 
M frequency values between 0 and f c . The number of segments in the data set, K 
is defined as N/4M. It is desired that N and M be in powers of 2. 

The statistical properties are next computed. The parameters used are the 
mean, jx, standard deviation, a ( root mean square when the mean is close to zero), 
skewness, s, and kurtosis, k [139]. The mean for a data set, Cj(t), of length N, 
given by 


is 


1 N 

" = N^ Ci 


(4.3.6) 


Similarly, the moments are as follows 


a = 
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N U 
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N U L 
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and 


1 N 

= i y 


j= i 
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Cj ~ /i 


- 3.0 


(4.3.7) 


(4.3.8) 


(4.3.9) 


The kurtosis moment ratio is usually 3 for an ideal gaussian process. Thus for a 
normal distribution the value of kurtosis is zero. A negative kurtosis implies a wide 
distribution shaped like a loaf of bread and a positive value relates to a sharply 
peaked distribution. A broad sense of the degree of nonlinearity in the response 
data can be obtained by observing how far this value is from zero. 
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The mean is for the most part close to zero, although it may be significant in 
the cases where thermal and acoustic loading are combined. The value of <r from 
the above equation is checked against the root mean square value obtained from 
the power spectrum estimate. They should be close. 

The peak probability distributions for the responses will be presented in a 
few cases. The probability densities of all the positive peaks in the response time- 
histories are computed. A peak is treated as a value of magnitude greater than its 
preceding and succeeding values. This may not be an accurate measure of the peak 
probability, say, when snap-through motion occurs. The number of sample points 
per ‘cycle’ in the time-history must not be too sparse. 

Apart from these, a few figures showing the probability distribution functions, 
PDF, are presented. The PDF shows the density of distribution of random samples 
about their mean. Again, this is a graphical check of the skewness - the deviation 
from the classic ‘bell’ shaped distribution - in case of nonlinear responses. A plot 
of the load PDF will be displayed for reference. 
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Chapter 5 


RESULTS AND DISCUSSION 


Numerical results are presented for numerous cases to illustrate the problem 
analysis and formulation discussed in earlier chapters. The random load model 
to be used is first described. The accuracy of the present formulation is verified 
by comparisons against existing solutions. Mesh and modal convergence studies 
are first presented for an unsymmetric cross-ply laminate. Most of the numerical 
examples that follow are for a symmetric laminate. The modal analysis for this 
laminate is detailed before random response results are presented. A modal conver- 
gence study under random load is carried out for this baseline laminate. Variation 
of responses under increasing load levels is tabulated. Temperature distributions 
considered include a simple uniform temperature assumption and a complete gener- 
alized non-constant distribution. Responses under uniform normal loading, as well 
as grazing-wave loading are then compared. Additional data to display the effects 
of other parameters, such as, plate thickness, initial stress, etc, are presented. 

5.1 Random Load 

The random pressure time-history generated has a frequency response function as 
given by Eq.(4.1.2). In this study, the sampling interval is set at At = 10 s and, 
thus, f c = 5000 hz. The Nyquist frequency, / c , easily includes the lower natural 
frequencies. Furthermore, the sampling interval, or the time-step size, is small 
enough to maintain low amplitude decay and period elongation. The number of 
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time-steps or samples obtained is set to N = 2 15 — 32768. The pressure spectrum 
has a bandwidth up to /„ which is chosen to be 4900 hz, and // = l/NAt = 
0.30517 hz , is the lower cut-off frequency. The value of f u was set at 4900hz to cover 
the whole r an ge of responses possible. A single filter was found to be sufficient. The 
spectrum of the resulting pressure time-history at 130db sound pressure level (SPL) 
is given in Fig. 4. Its probability density histogram is shown in Fig. 5. 

5.2 Verification of Formulation 

The finite element formulation presented in this study is first evaluated for accuracy 
using previously published results. The various intermediate stages in developing 
the program, were checked using simple free vibration and static nonlinear bending 
results. Linear forced vibration solutions (deterministic) were compared against 
solutions obtained from exact formulas. This helped verify the numerical integration 
algorithm. The final verification step is the only one detailed here, involving linear 
and nonlinear random vibration responses. 

5.2.1 Normal Random Load at AT = 0 

Comparisons of results from the present formulation are made against those from 
Ref. [89] and Ref. [105]. Mei and Prasad, in Ref.[89], compute, using an analytical 
Galerkin single-mode approach and the equivalent linearization technique, random 
response solutions for square simply-supported symmetric cross-ply laminates of 
various length- to- thickness (a/h) ratios. Results with and without transverse shear 
effects have been listed. The classical formulation of the governing equations uses 
the Airy stress function, thereby neglecting membrane inertia terms. All com- 
parisons made here are for transverse shear deformable behavior. The material 
properties for the graphite-epoxy laminate used here are as follows 
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Fig. 5 PDF of Acoustic Pressure 



E 2 = 0.75 x 10 6 psi 
E x fE 2 = 40 
G 23 /E /2 ~ 0-5 
G 12 /E 2 — G 13 / E 2 = 0.6 
v \2 = 0.25 

p = 2.4 x 10 -4 lb - s 2 /in. 4 
a = 12m., a/h = 200 

damping ratio, £ = 1%, with modal damping, Cj 2<^o;,. 

A shear correction factor of 1C = 1 is used, to coincide with. the unit ‘tracing 
constant’ used in Ref. [89]. The boundary conditions used are for a plate simply- 
supported on all sides, with immovable inplane boundary conditions as shown in 
Fig. 6. Since the plate is symmetric, consisting of cross-ply layers only, a 4 x 4 
quarter-plate model is used. The analysis leading to the mesh size used here is 
detailed in section 5.3.1. The fundamental frequency obtained is at 69.992 hz, which 
exactly matches the frequency corresponding to the nondimensionalized parameter 
u 0 = 18.88 given in Ref. [89]. 

The nonlinear simulation to be done requires the use of membrane mode-shapes. 
Of the lowest few membrane modes, the first mode was found to be inactive, i.e., 
its presence made no difference to the deflection rms values. Therefore, the next 
3 membrane modes are used. As seen in results below, sufficient agreement in the 
deflection results is obtained with these three membrane modes. Thus, simulation 
is performed with a total number of modes R = 4, with 1 bending mode and 3 
membrane modes. 

Numerical values are compared at 130db SPL, for < 7 , a nondimensional root- 
mean-square (rms) deflection parameter for the maximum deflection at the center 
of the plate, given by 
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Fig. 6 Simply- Su] 



_ _ WcrErf* 

a 4 y^PSD fi 

where a is used to represent the standard deviation for maximum deflection, (same 
as rms since the mean is negligible), f\ denotes the fundamental frequency, and 
PSD stands for the magnitude of the power spectral density ( psi 2 jhz ) for 130db 
SPL (see [89]). Using a reference pressure level value 
Pre/ = 2 x 10 -5 Pa = 2.9 x 10 ~ 9 psi 

we have, 

PSD = p 2 Tt j^ SPL ^ Q — 8.425 x 10 -5 psi 2 /hz 
The results for linear and nonlinear responses are given in Table 1. A graphical 
representation of responses at various load levels is shown in Fig. 7. As shown in 
the table, the difference from Ref. [89] is within 4%. 



linear 


nonlinear 



Ref. [89] 

present 

Ref. [89] 

present 

a (in.) 

0.3955 

0.3800 

0.0955 

0.0917 

a 

0.4024 

0.3860 

0.0972 

0.0933 

a jh 

6.592 

6.333 

1.592 

1.53 

Difference 

- 

3.8% 

- 

3.7% 


Table 1. Cross-Ply Laminate Rms Deflection 

There seems to be very good agreement between the results in both Table 1 as 
well as Fig. 7. Some of the factors in the computation procedure that differ between 
the two methods need to be highlighted. The Lagrangian quadratic element used 
in this study, is known to converge to the exact deflection from below (frequency 
from above). The numerical integration algorithm as well as the Newton-Raphson 
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iteration technique are known to produce conservative results. However, this influ- 
ence is rather small for a system with moderate nonlinearity (note that ojh « 1.6). 
These are, however, factors that could predict lower values than a classical equiva- 
lent linearization method. A point to mention is the level of membrane contribution 
for this case. The ratio of membrane modal rms values to the fundamental mode 
rms value were in the order of 10 . 

At this stage, the numerical simulation formulation has been verified for non- 
linear deflection responses. This confirms the accuracy of the formulation in the 
absence of thermal loads. A numerical comparison of strain data is presented next. 

A verification for the strain values was made against data from Ref.[105]. In 
Locke’s study, a thermally buckled isotropic plate under a symmetric acoustic load 
is analyzed using a post-buckling formulation in conjunction with the equivalent 
linearization method. The membrane stiffnesses were eliminated by static reduc- 
tion. A quarter-plate finite element model with 36 24-degree-of-freedom rectangular 
elements was used. The results for simply-supported boundary conditions with im- 
movable edges, as in Fig.6, are used here. The properties of the plate are (with a 
shear correction factor included for the present formulation) 

E = 10.5 x 10 6 psi 
v = 0.3 

p = 2.588 x 10 -3 lb — s 2 /in.* 
a x b x h = 15xl2x 0.04 (in.) 

-C = 1% 

K = 7T 2 / 12 

/! = 44.078 hz 

Single-mode finite element results were compared with classical solutions by 
Locke in Ref. [105]. Verification here is for a plate with no thermal load. The 
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finite element model used here again has a quarter-plate 4x4 mesh. The numerical 
simulation used the lowest few membrane modes, of which the first mode is relatively 
weak; its influence on the deflection results being negligible. The next 2 membrane 
modes were used along with one bending mode, i.e., 

R = 3 ; 1 bending mode and 2 membrane modes. 

The single-mode nonlinear random responses at 95db and 125db SPLs are 
compared below in Table 2. The maximum strain is in the y-direction and its rms 
value is denoted by £ rmJ , given in micro-strain units, 10 -6 in. /in. The rms and 
standard deviation for strain Eire not identical, since the mean is no longer small. 



95db 


125db 



ajh 

^rms (/^) 

a/h 

^rms (/^ e ) 

Classical [105] 

0.391 

22.41 

2.771 

362.95 

FEM [105] 

0.390 

22.58 

2.766 

366.49 

present 

0.380 

21.42 

2.6925 

301.82 

Difference 

2.77% 

4.4% 

2.8% 

16.8% 


Table 2. Rms Deflections and Strains for Isotropic Plate 

The difference above has been computed against the equivalent linearization 
classical solutions given in the first row of the table. The rms deflection values 
agree very well and so does the strain at 95db. The difference in strain values at 
the higher load has possible cause in a few factors. The differences in the membrane 
formulations between the two methods is a possible factor, since membrane action is 
stronger at high load levels. The interpolation functions and the finite element mesh 
used also influence the displacement derivatives. The displacement and rotation 
interpolations in this study are C° continuous, resulting in their derivatives being 
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discontinuous across element boundaries. In Ref.[105], the bending displacement 
function has C 1 continuity. The strain computations are, therefore, likely to be more 
accurate. The inplane displacement functions used are C° continuous. But the same 
interpolations are again used to calculate the strains and curvatures, making these 
continuous. The figures posted here for strains are actually the maximum elemental 
strain values, at the centroid of the element and not at a particular node. In effect, 
an average over an element is being compared with the maximum occurring at one 
of its nodes. 

5.2.2 Normal Random Load with Constant AT 

The verification for this case is done against the thermal post-buckling solutions 
of Ref. [105]. Equivalent linearization results for single- mode rms deflection and 
strain responses are generated at multiples of the critical buckling temperatures. 
The present simulation results are compared against these, at the critical buckling 
temperature, A T cr . 

The lowest non-zero eigenvalue of the buckling eigenvalue problem corresponds 
to A T cr , scaled by the AT value used to build up the thermal stiffness matrix, 
A'at (see Eq.4.2.6). The value of A T cr , using the same material properties and 
boundary conditions listed in section 5.2.1, was found to be A T cr = 0.9° F. Again, 
a single bending mode is used with two membrane modes in the simulation. Using 
a uniform temperature distribution at AT = A T cr and a 95db acoustic load, the 
rms values obtained are listed in Table 3 below. Comparisons were made at a lower 
decibel load, since the effect of thermal loading is more pronounced than at higher 
levels. 

There seems to be reasonably good agreement given the differences between the 
two methods of formulation discussed earlier. A major difference between the two 
formulations that affect combined loading results, is in the post-buckling solution 
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used by Ref. [105] for both the classical and finite element methods. The static 
deflection shape obtained from the thermal loading problem is used as an initial 
deflection in the random vibration problem. The membrane action is accounted for 
by this static mode, making the model stiffer than a multiple-mode-assumption. In 
the case of the simulation, the thermal and acoustic loading occur simultaneously, 
which amounts to the thermal load also being treated as part of a dynamic problem. 
This is an important difference as it results in higher buckling eigenvalues being 
involved in the response. Their effect will be studied later. 



ajh 

^rms (A^) 

Classical Ref.[105] 

0.497 

29.22 

FEM Ref. [105] 

0.496 

29.42 

Present 

0.4965 

26.0 

difference 

0.1% 

11% 


Table 3. Rms Deflections and Strains for Combined Load 
5.3 Numerical Examples 

The following sections present responses for various acoustic and thermal load cases 
for a specified baseline plate model. A mesh and modal convergence analysis is first 
performed on an unsymmetric laminate. The baseline plate used is a symmetric 
laminate. Results axe presented for response under varying decibel levels. Fur- 
ther simulation for various temperature distributions, the grazing incidence acoustic 
load, etc., is performed at a baseline load level of 130db SPL. 


69 



5.3.1 Mesh Convergence 

Different mesh sizes are tested to represent a full-plate model of an unsymmetric 
laminate. The changes in the natural frequencies, nonlinear random deflections and 
strains are observed. An unsymmetric laminate is used, since separate membrane 
modes are not necessary when Bij 7^ 0 (see section 4.2.2). The unsymmetric lam- 
inate used is a two-layer cross-ply graphite-epoxy panel of a (0/90) lay-up. The 
plate is assumed fully clamped on all edges with immovable inplane displacements, 

i.e., 

U = V = w = xl>x = V’y = 0, on all 4 edges 
The material properties for the laminate are 

Ei = 22.5 x 10 6 psi 
E 2 = 1.17xl0 6 psi 
G 2 3 = 0.66 x 10 s psi 
Gi 2 = G 13 — 0.4 x 10 s psi 
v \ 2 = 0.22 

p = 1.45S x 10 ~*lb-s 2 /in. 4 
a x b = 15 x 12 (in.) 
h = 2 x 0.024 (in.) = 0.048 (in.) 

c = 2%, c, = 

K = tt 2 /12 

Free-vibration eigenvalues and eigenvectors are obtained for three different 
mesh sizes. The mesh sizes used are a 6 x 6, 8 x 8 and 10 x 10 grid of rectangular 
45-degree-of- freedom elements over the full plate. The number of mode-shapes re- 
quired to obtain converged nonlinear rms maximum deflection and rms maximum 
strain values were determined using the 8x8 mesh model. The details of this modal 
convergence study will follow in section 5.3.2. A total of seven modes ( R = 7) is 
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used for nonlinear random simulation at 130db SPL for all the three mesh models. 
The highest of these modes is mode (4,1). Table 4 below shows results for three 
parameters, namely, the frequency for mode (4,1), the center deflection, <r, and the 
maximum strain (y-axis), e rmj . The strain values are, as mentioned before, given 
at the element centroid. They represent am averaged value over the area of the ele- 
ment. With a changing mesh, the location of the element where strain is maximum 
changes. Results from the 10 x 10 mesh are used as a basis for comparison of results 
from the smaller meshes. This mesh was selected as a reference mesh, based on 
linear random response convergence studies by Robinson [143]. 


mesh size 

mode (4,1) 

(7 

trms 

6x6 

4.4% 

3.8% 

14.3% 

8x8 

0.7% 

-0.8% 

-10.2% 


Table 4. Mesh Convergence for (0/90) Plate 

The minus signs above indicate that the values were lower than the 10 x 10 mesh 
values. As can be seen, mesh convergence for the nonlinear quantities, deflection 
and strain, is not monotonic. Robinson [143] shows monotonic mesh convergence in 
linear random deflection response. Locke [144] conducted mesh convergence studies 
for thermal post-buckling nonlinear deflections and stresses for a clamped isotropic 
plate. All these converged monotonically. The added change in the present case is 
the modal transformation combined with nonlinear random simulation. Therefore, 
the possible factor in the present convergence behavior is the modal cross- correlation 
in the transformed equations of motion. 

The 6x6 mesh shows a fair amount of difference even in the natural frequency 
results. Therefore, this mesh will not accurately model the required mode-shapes. 



The 8x8 mesh shows excellent comparison with regard to both the frequency and 
the nonlinear deflection. This indicates that there is adequate resolution in the 
modelling of the required mode-shapes. The strain value shows poorer convergence 
behavior, but the difference is still in the 10% range. Amongst the various stresses 
due to thermal loading studied by Locke [144], the membrane stress at the edge 
midpoint showed the poorest convergence behavior. Given adequate resolution in 
modal representation, the strain values can be expected to vary with element sizes, 
as they represent averages over progressively smaller areas. 

Th 10 x 10 mesh size, resulting in 1805 degrees-of-freedom, is a large problem 
for the eigenvector computation method that was used here. In this regard it is 
preferable to use the 8 x 8 mesh. The difference in strain data between the two is 
within a reasonable range. Hence, for the most part, the 8x8 mesh is chosen in 
the following examples. 

5.3.2 Modal Analysis 

The modal analysis leading to the random results in the above section is detailed 
here. The same analysis for the symmetric laminate to be used is also presented 
in this section. This symmetric laminate has no coupling between bending and 
membrane mode-shapes and, therefore, will involve modal convergence studies for 
both the membrane and bending mode-shapes. 

Modal convergence for the unsymmetric (0/90) laminate is carried out using the 
8x8 mesh on the full clamped plate. The lowest ten mode-shapes in ascending order 
of their natural frequencies, axe as follows: (1,1), (2,1), (1,2), (2,2), (3,1), (3,2), (1,3), 
(2,3), (4,1), and (4,2). Their frequencies ranged from 75.5 hz for the fundamental 
frequency to 421.6 hz for the tenth natural frequency. Random response simulations 
are run at 130db SPL, with increasing number of modes included. The modal 
rms values are computed along with rms maximum deflection, at the center, and 
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maximum strain, along y-axis. The element corresponding to the maximum strain 
is, for this mesh and load, not at the long edge midpoint. It was found to be one 
element away from the edge. This can be expected, since the maximum strain for 
a clamped plate is ‘close’ to the edge. A refined mesh picks this maximum value 
at an element away from the edge. Further discussion on this follows under section 
5.3.3. 

Figure 8 shows the modal convergence curves for deflection and strain for this 
case. Of all the lower modes listed above, no response was detected at modes (1,2) 
and (2,2), the third and fourth modes, respectively. The modal time-histories at 
these two modes showed rms values below 10 — 12 . Further, including or excluding 
them did not affect the deflection and strain responses. The deflection as well as 
strain data show reasonable convergence within the first four active modes itself. A 
more confident estimate of converged values is with seven modes, as there is modal 
activity, albeit weak, at modes beyond these. 

The converged rms maximum deflection and strain values for this case are 
a/h = 1.737 and e rms = 377.3fie. 

Baseline Plate 

The modal convergence studies for the baseline plate are discussed next. The prop- 
erties of the plate are the same as those listed in section 5.3.1. Its symmetric lay-up 
consists of eight layers, each 0.006m. thick, in a (0/ + 45/ - 45/90), configuration. 
All other parameters, including the total thickness, remain the same. The plate is 
assumed fully clamped. Two finite element meshes are generated for this plate, the 
8x8 mesh and a smaller 8x4 mesh. This smaller mesh has been chosen with the 
grazing load model in mind. Since the pressure-wave is assumed uniform along the 
y-axis, the number of elements along this direction is reduced. In all the results to 


73 




74 


Fig. 8 Modal Convergence for 0°/90° Clamped Pi 


follow, the mesh and modal model used is specified for each case. 

Free vibration eigenvalue problems were solved for both these meshes. An 
inspection of the bending modes indicated the presence of some ‘skewness’ in their 
shapes, despite the symmetric lay-up of the laminate. For instance, the fundamental 
(1,1) mode was found to be not perfectly symmetric. Neither are the anti-symmetric 
modes perfectly anti-symmetric. The (2,2) mode deviates from the classic single sine 
function curve as we approach the centerline of the plate. This can be observed in 
the figures for the first four modes plotted in Figs.9 and 10. An observation of 
the centerline for these mode-shapes [except mode (1,1)] shows the deviation from 
expected conventional shapes seen in isotropic or orthotropic laminates. 

This type of skewness should be expected since the baseline laminate has angle- 
ply layers at ±45°. These angle-ply layers result in non-zero values for D l6 , D 2 6 
and their transpose positions (the twist-bending coupling terms) in the laminate 
bending stiffness matrix, D{ } (Eqs.2.3.9-11). Although the Ai6, A 2 6 positions of 
the extensional stiffness matrix are also non-zero, they are maximum for a three- 
layer angle-ply laminate and drop-off sharply as the number of layers increase. In 
keeping with this property, the A ]6 , A 26 values were found to be negligible. A 
check of the membrane mode-shapes showed an absence of skewnesss, as there is no 
extension-twist coupling. The drop-off in the stiffness values with increasing number 
of layers is not so steep for the twist-bending coupling terms. They vary as L 2 /L 3 
for tin L-layered laminate. The influence of this twist coupling on the skewness of 
the mode-shapes can be illustrated by redoing the eigenvalue analysis for the same 
laminate with D i6, D 2 6 and its transposes set to zero. Fig. 10 shows the (2,2) mode 
with and without twist coupling. The skewness has evidently vanished. Due to the 
presence of this skewness, all the lowest bending modes [(1,1), (1,2), (2,1), (2,2), 
(1,3), (3,1), . . .] up until the twelfth natural frequency are saved for convergence 
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Fig. 10 Mode (2,2) With and Without Twist-Coupling 


77 


studies. As the mode-shapes get more complicated [(3,3), (2,4), etc.], it is difficult 
to accurately label them so. Strictly speaking, it is rather erroneous to assign, say 
the mode (3,3), as done here, since the skewness seriously distorts its shape over 
the plate. The membrane mode-shapes are combinations of similar shapes for the 
x and y-direction displacements, u and v, respectively. For example, the lowest 
mode to be used here has shape (1,2) in x-direction and (2,1) in the y-direction. 
Convergence of nonlinear random deflection and strain responses are now done. 
In the interest of accuracy, all data presented are for the 8x8 elements mesh. 
Bending as well as membrane modes need to be picked. The modal convergence 
for number of membrane modes required is done using one bending mode, namely 
the first mode (1,1). A total of twelve membrane mode-shapes were saved ranging 
in frequency from 10207hz to 22336 hz. (It must be clarified here that, it is not 
the response solutions at these frequencies that are of interest, but the effect of 
these modes on the bending mode responses. This effect is represented in the 
modal equations through the nonlinear bending-membrane coupling tensor, kl pra 
(Eqs.2.5.13). Furthermore, the level of modal response at these frequencies is far 
smaller than those at the bending mode frequencies. This is why the excitation 
spectrum is not extended up to these high frequencies.) The highest mode amongst 
these has shapes (2,3) and (1,2) in the u and v displacements, respectively. Of 
these twelve membrane modes, six were found to be inactive, i.e., their modal rms 
values were zero and the deflection and strain results remained unaffected by their 
exclusion. Convergence results are with increasing number of active modes included. 

Values for rms maximum deflection per unit thickness (cr/h), and maximum 
strain (e rm ,), at 130db SPL for increasing number of membrane modes are plotted 
in Fig. 11. The maximum strain was found to be along the y-direction. The strain 
values are specified in micro-strain units. With the 8x8 mesh chosen, the element 
showing maximum y-axis strain at this load level is not at the long-edge mid-point, 
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Fig. 11 Membrane Mode Convergence for Baseline Plate 


but at the next element towards the plate center, one row away from the edge. The 
sixth mode in the figure, where solutions have converged, in fact corresponds to the 
twelfth membrane natural frequency. 

Bending mode convergence is studied next using linear and nonlinear simu- 
lations at 130db SPL. Nonlinear simulations use six membrane modes from the 
convergence shown above. Results for linear and nonlinear center deflection and 
y-axis strain are presented in Fig. 12. As in the case of membrane modes, not all the 
lowest bending modes are active. The mode-shapes that were found to be active 
correspond to (1,1), (2,2), (1,3), (3,1), (3,3) and (2,4). Their natural frequencies 
are 102.3, 300.12, 338.45, 417.4, 594.0 and 652.9 hz, respectively. In the case of the 
linear simulation, as seen in the figure, even single-mode solutions are fairly close 
to converged four-mode results. The convergence in the nonlinear simulation is not 
monotonic for both the deflection and strain curves. This is due to the presence of 
coupling between modes. This does indicate that the practice of neglecting cross- 
correlation terms in the modal equations, as done in most earlier analytical work, 
may lead to erroneous inferences on plate responses. 

The converged rms values for a jh and e rni , are 6.4 and 628.0 j/e, respectively 
in the linear case, and 1.86 and 438.5 fie, respectively for the nonlinear simulation. 

Thus, the nonlinear solutions for the baseline plate at 130db SPL can be said to 
have converged with six membrane and six bending modes ( R = 12). The various 
cases analyzed in the succeeding sections, use this set of modes for the 8x8 mesh. 
Since the center deflection value shows better convergence characteristics both in 
mesh size and number of modes than strain, the smaller 8x4 mesh with fewer modes 
has been used where only deflection data axe presented. This smaller baseline plate 
model uses just two membrane modes based on the membrane convergence in Fig. 11 
and four bending modes based on the bending convergence in Fig. 12. 
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nee for Baseline Plate 


5.3.3 Load Variation 

The next set of results displays the changing qualitative characteristics of panel 
responses under increasing sound pressure levels. The decibel levels used ranged 
from 90db SPL, where the deflections can be expected to remain within the linear 
regime of panel response, to 130db SPL, where as seen previously large deflections 
have set in ( a/h w 2). All the load levels had spectra similar to that in Fig. 4. 
The pressure rms value obtained by integrating the frequency response function is 
simply scaled to suit the required SPL. The simulations were carried out with the 
8x8 mesh plate model. 

The variation of deflection response with SPL is shown in Fig. 13. The hard- 
ening type nonlinearity is seen in the reduced nonlinear deflection values at higher 
loads. Data using the smaller baseline mesh defined earlier is also plotted in this fig- 
ure. This gives a good estimate of the accuracy of a coarser mesh and fewer modes 
(R — 6), in predicting nonlinear deflection values. As can be seen the smaller mesh 
is fairly adequate as far as deflection values concerned. All other data in this section 
pertain to the 8x8 mesh model. 

The effect of increasing load, and thereby, increasing nonlinearity, can also be 
seen in the power spectra for the center deflection shown in Fig. 14. The spectrum 
for 90db SPL load shows distinct peaks corresponding to the (1,1), (3,1) (1,3) and 
(3,3) modes, respectively. The mode (2,2) which does show up in the modal rms 
values, here is a very small peak. The mode-shape in Fig. 10 for this mode has a very 
small magnitude at the plate center. The spectra at loads HOdb abd 130db show the 
hardening type nonlinearity, in the shift in their response peaks. There is significant 
broadening in these peaks, another indicator of nonlinear random response. 

Table 5 shows the various statistics for rms maximum deflection at varying 
load levels. The standard deviation and rms coincide since mean deflection is zero. 
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Fig . 



The skewness and kurtosis, denoted by s and k, are non-dimensional quantities 
listed in columns 3 and 4. 


SPL (db) 

ajh 

s 

k 

90 

0.0638 

0.0558 

0.614 

100 

0.1969 

-0.004 

0.32 

no 

0.524 

-0.0026 

-0.743 

120 

1.147 

-0.008 

-0.668 

130 

1.86 

0.0152 

-0.796 


Table 5. Maximum Deflection Statistics vs Sound Pressure Level 

For practical purposes the skewness values are very small and indicate a rea- 
sonably symmetric probability density distribution about the mean. Graphical il- 
lustrations of the probability density function (PDF) for the nonlinear deflection 
response at 90db and 130db load levels are displayed in Figs. 15 and 16. The PDF 
of the linear response at 130db is shown in Fig. 17 for comparison. As can be seen 
from these figures, the PDF is pretty symmetric about zero. A good indication of 
the increasingly high density of large magnitude responses is the trend in the kur- 
tosis values. As mentioned in Chapter 4, a negative kurtosis indicates a bread-loaf 
shape of the PDF and a positive kurtosis implies a narrower distribution peaking 
towards the center. The kurtosis numbers in Table 5 do reiterate the behavior seen 
in the figures. This nonlinearity as indicated by the kurtosis value at 130db is again 
reaffirmed in the peak probability density function for the large deflection response 
shown in Fig. 18. Only positive peaks in the deflection time-history have been con- 
sidered. A plot of the same for the linear deflection at 130db SPL is also shown in 
Fig. 19. A comparison of the two shows a marked shift in the peaks towards higher 
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Fig. 16 PDF of Nonlii 
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Fig. 18 Peak PDF of Nonlinear Deflection at 130db SPL 
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magnitudes in the nonlinear case. The Rayleigh distribution of the linear response 
peak density has also widened considerably. 

r J The behavior of the normal strains in the plate with increasing load levels 
is next analyzed in a similar fashion. However, in this case, a breakdown of the 
various components in the maximum strain is presented. This clearly displays the 
changing quantities in the total strain with large deflection behavior setting in. 
The maximum strain was found to be in the y-direction (the shorter dimension), 
with x-direction strain reasonably well-behaved throughout. Therefore, the focus 
is on the maximum y-direction strain, Cy. The maximum is computed at thickness, 
z — h/2. Unlike in the case of maximum deflection, which is at the center of the 
plate, the location of the maximum strain shifts with increasing load, for the 8x8 
finite element mesh considered. The maximum strain is close to the mid-side of the 
long edge for all the load cases, since this is a clamped plate. 

The power spectra at different loads for this strain are displayed in Fig.20. 
As in the case of deflection, the shift as well as widening in the response peaks is 
present here, too. The strain response at 130db SPL has almost reached a uniform 
broadband behavior, across the frequency spectrum. The power spectrum for strain 
at 90db generated using the smaller 8x4 mesh with a total of six modes, as 
mentioned before, has been included to draw attention to one important observation. 
A comparison of the 90db strain response spectra for this mesh with that for the 
8x8 data shows a difference in the peak magnitudes. The higher magnitude for 
the second peak in the smaller mesh data indicates a greater contribution to strain 
at this mode. This indicates that this coarse mesh is too stiff to accurately predict 
strain contribution from this (1,3) mode. 

Table 6 lists rms value of all the components of rms y-strain, e rm3 - e m , the 
linear membrane normal strain vector; £*, the linear strain vector corresponding 
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Fig. 20 Power Spectrum of NL Y-strain at Different SPL 



to curvature and e w , the nonlinear strain vector due to large deflection (see Eqs. 
2.2.6-10). The mean value of e denoted by e M , is also included. All units are in 
micro-strains, /xc. The skewness and kurtosis values axe not presented for strain. 
They follow the trend in the deflection statistics but tend to be more extreme. All 
the statistics for the curvature component were found to be indicative of a zero-mean 
gaussian process at all loads. 

At lower decibel levels, the dominant component in the total strain is the 
linear strain due to curvature, c K . This is maximum at the edge. As the load 
level increases, the nonlinear strain component, e w , begins to gain magnitude. This 
component is highest a little away from the edge of the plate. This can be expected, 
since, right at the clamped edge, membrane extension is constrained. Therefore, at 
120db and 130db levels, the element for which the strain reached a maximum has 


SPL 


^ rma 

m 

rmj 

e K 

rmj 

e w 

rmj 

90 

0.2285 

6.248 

0.035 

6.245 

0.4324 

100 

2.21 

19.54 

0.477 

19.09 

4.35 

110 

17.65 

55.53 

3.42 

49.21 

27.91 

120 

97.6 

182.4 

4.5 

105.6 

151.7 

130 

254.8 

438.5 

15.8 

217.95 

389.4 


* all in micro-strains 

Table 6. Rms Strain Components vs Sound Pressure Level 


moved one row towards the interior of the plate. 

The PDF distributions for this y-axis strain at 90db and 130db load levels are 
shown in Figs. 21 and 22. The skewness in the strain at 130db is clearly noticeable. 
The peak PDF (positive peaks only) for 130db nonlinear strain response is given in 
Fig. 23. Again, for comparison the linear peak PDF is displayed in Fig. 24. 
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Fig. 21 PDF of Nonlinear Y-strain at 90db 
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Fig. 22 PDF of Linear and NL Y-strain at 
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Fig. 23 Peak PDF o£ Nonlinear Y-strain at 130db SPL 
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Fig. 24 Peak PDF of Linear Y-strain at 130db SPL 


As can be observed in the PDF figures, the increase in skewness is not accompanied 
by a strong deviation from a gaussian distribution. There is a slight move towards 
a Rayleigh distribution type shape in the nonlinear strain PDF. The wider distri- 
bution of strain in the nonlinear PDF is reflected in the peak density distribution 
of Fig. 23. The broader range of peak magnitudes in the nonlinear peak PDF is 
concentrated in a smaller range in the linear response peak distribution of Fig.24. 
Further, there is a shift towards higher magnitudes in the nonlinear peak distri- 
bution. This shift is not as pronounced as in the case of deflection (Fig. 18). The 
‘random’ characteristics of nonlinear strain response are, for the most part still de- 
termined by its linearly-behaving curvature component. The nonlinear membrane 
action contributes to the strain primarily as a mean value increase. 

5.3.4 Variation of Temperature 

The changing responses under increasing temperature loads are now analyzed. The 
baseline plate is used with load at 130db SPL. The smaller 8x4 mesh is used here, 
since the results presented are for deflection behavior only. The total number of 
modes used to approximate the deflection is R = 6, as described before in section 
5.3.2. The comparison in Fig. 13 has shown that the center deflection rms values are 
reasonably accurate even with this model. The lowest few buckling temperatures 
were determined for this baseline plate configuration using Eq.4.2.6. The values 
obtained are 

A T er , = 37 °F 
AT C r 2 = 49 °F 
A T crs = 83 °F 

Simulation results are generated at increasing levels of AT, ranging from 0° 
to 80°F. A toted of 15 different temperatures are considered, keeping in mind the 
buckling values above. These results combine both the dynamic random response as 
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well as the thermal response. The variation of rms deflection (<r) with temperature 
is shown in Fig. 25. Strain results, although not presented here, did follow the same 
pattern of behavior. The peaks in the deflection response, denoted as ATi, AT-* and 
AT3, approximately correspond to the above three buckling temperatures. When 
the same simulations were done with a single bending mode, the peaks more closely 
matched the buckling temperatures. These peaks are indicators of reduced modal 
(vibration modes) stiffnesses, at those temperatures. As the temperature passes 
the buckling value, the response increases and then decreases as the panel stiffness 
starts increasing again. These peaks are results for a system that includes the 
nonlinear stiffnesses, the inertia and damping terms. Therefore, it can be expected 
that the temperature values at which they occur do not exactly match the buckling 
eigenvalues given above. 

This multiple peak characteristic of response under combined loading is not 
identified with a thermal post-buckling analysis. The buckled plate is modelled by 
the equivalent of its critical buckling temperature mode-shape. This eliminates the 
singularities that occur in the [A" - K AT ] system at higher eigenvalues. 

Snap-Through Motion 

At temperatures much higher than those used above, snap-through motion in the 
large-amplitude deflection has been reported ([109], [125], [134]). A simulation 
at 130db SPL with temperature increased to AT = 350°F, did result in such a 
deflection response. The baseline plate model has six modes in an 8 x 4 mesh. 

The time-history for maximum deflection, w max (t), at snap-through is dis- 
played in Fig. 26. A very short segment of about 0.1 sec , out of a sample of over 
3 sec, has been extracted. The ‘oil-canning’ phenomenon, where vibration occurs 
about some non-zero amplitude for short bursts of time is clearly visible. The 
parallel segment for nonlinear time-history at AT = 0 °F, is shown in Fig. 27 for a 
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Fig. 25 Nonlinear Deflection vs uniform Temperature at 130db 
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comparison. A linear simulation time-history at the same load and with no temper- 
ature is shown in Fig. 28. This figure shows a much slower moving curve since the 
response is dominated by low frequency oscillations. A comparison of the former 
two with Fig. 28, gives an idea of the large high frequency content in the nonlinear 
response. 

The power spectrum for the snap- through response is, in fact, very similar to 
the 130db response spectrum with no thermal load (Fig. 14), except for some very 
low frequency content. This is due to the snap-through oscillation being picked up 
as a very slow-moving mode. 

Uniform Temperature at Different Loads 

The response data for nonlinear deflection at 130db load under uniform tempera- 
ture, does not deviate much from the zero temperature response values at the same 
random load. Earlier research has stated that thermal loading effects are dominant 
at lower random load levels [83]. This aspect of plate behavior is now studied. 

Two load levels are used - 130db and 90db SPL. Random response time- 
histories at these loads with AT = 80°F are generated. For the sake of better 
accuracy in strains, the larger model of twelve modes ( R = 12) in the 8x8 mesh is 
used to generate data. Table 7 below lists rms deflection and strain values for both 
these load cases, at AT = 0° and 80° F. The statistics include the non-dimensional 
mean and rms deflection values and the mean and rms values for the maximum y- 
axis strain, and its curvature component, e K . This component has been included 
to show its reaction to thermal load. 

The maximum strain location is different for these two load cases. In the case 
of 130db, it is in the element close to mid-point of the long-edge one row away from 
the edge. In the 90db case, maximum is at the element on the edge itself. 
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130db 

90db 



0°F 

80°F 

O^F 

80° F 

H Wmax/h 

-0.007 

-0.01 

-4.31 x 10“ 5 

-1.101 

rms w max /h 

1.86 

1.94 

0.0637 

1.103 


254.8 

278.1 

0.228 

165.07 

^rms 

438.47 

479.12 

6.25 

165.6 


-0.527 

-1.3 

0.0043 

111.9 

frm, (/*«) 

217.95 

244.256 

6.245 

112.05 


Table 7. Statistics at 90db and 130db SPL 

The influence of the thermal load at 90db is very pronounced, whereas at 130db 
there is a very small difference. Both the deflection and strain show a high mean 
at 80°F. A big increase in the linear component, e K is also observed. 

The change in the responses at the lower load level is also shown in Figs. 29 
and 30. The first is a comparison between the power spectrum for the nonlinear 
deflection with and without thermal load at 90db SPL. The spectrum in this figure 
can be understood by the time-histories plotted in Fig.30. The vibration at 80°F 
takes place about some buckled position. This type of change is not present at 

130db SPL. 

Non-uniform Temperature AT = T(x,y,z) 

The thermal load modifications to model a non-uniform temperature distribution 
axe now implemented. Data from experimental setups suggests that a panel exposed 
to a bank of heat lamps intended to create a temperature distribution as close to 
being uniform over the plate surface as possible, is, in reality, not at a uniform 
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temperature. There is a drop in temperature near the edges of the clamped panel 
[145]. This type of distribution, which is symmetric over the z -y plane, is modelled 
here. Figure 31 shows the temperature variation assumed along both the x-axis and 
the y-axis for the baseline plate. 

The gradient through the thickness of the plate is known to be quite small. 
Since the plates are thin, and the duration of exposure is sufficiently long, this 
seems reasonable. This gradient usually does not exceed about 20 F. Hence, a 
linear variation along the z-axis is assumed, with the maximum being on the external 
surface (Fig.3). Therefore, the expression for the temperature variation in the plate 


can be written as 

A T{x,y,z) 


T(x) T(y) [ T a + 



where T 0 is a constant, corresponding to the temperature increase above ambient 
atmospheric temperature. The slope in the z-direction is given by T\/h. 

As mentioned in the introduction, this study is motivated by the structural 
design requirements for the proposed high speed civil transport aircraft. It is es- 
timated that the operating temperatures of interest are at about 350 F . Ambient 
temperatures are assumed to be 70 °F. Therefore, if the temperature differential for 
the outside surface of the panel is at 280° F, the mid-plane temperature differential 
is set at 270°F, to accomodate the gradient. Thus T 0 = 270°F and T\ = —20.0. 
The surface distribution, T(z,y) = T(x)T(y), is simply scaled by the T 0 value as- 
sumed. For a uniform temperature distribution, these are set to unity and Ti is set 


to zero. . 

Since the previous result of this section showed a small thermal load influence at 
the 130db baseline load level, this non-uniform temperature analysis is carried out 
at the lower 90db SPL load. The finite element model uses the baseline 8x8 mesh 
with twelve modes. The effect of each of the temperature distributions, namely, 
uniform temperature, AT = T 0 \ temperature with surface variation only, 
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Fig. 31 Components of Non-uniform Temperature Distribution 
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AT = T(x,y); and temperature with variation on the surface and through the 
thickness, AT = T{x,y,z ) is studied. Results for mean and rms deflection and 
strain are compared for these three themal load assumptions. Table 8 lists the rms 
values for this case. The element with maximum strain is one row away from the 
edge in this case. 

Overall, the differences in the temperature loadings do not seem to have affected 
the response characteristics to a great degree. A quantity to observe, though, is the 
curvature component of the strain, e K . It shows a significant increase from a uniform 
temperature distribution to a surface gradient distribution. (In fact, at the very 
edge of the plate, this component was found to be larger than the numbers above, 
although the total strain still remained lower.) The temperature distribution, T(y), 
shown in Fig.31, has the largest gradient along the y-direction close to this edge. 
This surface gradient can be expected to create high thermal strain in this vicinity. 
This possibly explains the change in e K seen above. 



To 

T(x,y) 

T{x,y,z) 

ft ^mai / ^ 

2.575 

2.53 

2.53 

rms W m ax/h 

2.58 

2.53 

2.58 


578.86 

581.7 

581.7 

(rmi (A^) 

580.6 

582.7 

582.7 

(/*«) 

76.6 

103.3 

103.3 

C, ) 

514.6 

489.1 

489.1 


Table 8. Statistics vs Temperature Distribution - 90db 


110 



fcS 


ca 


(£3 



m 


The inclusion of a gradient of this particular magnitude in the z-direction does 
not have any effect on the response statistics. The motion of the plate is so domi- 
nated by the thermal force, that the thermal moment countering the moment due 
to the dynamic load, makes very little difference. This may not be the case for 
higher gradients. 


5.3.5 Grazing Load 

In this section results are presented for the pressure modelled as a plane- wave trav- 
elling along the length of the plate. It is assumed that the angle of incidence of 
Fig. 3 is fixed at 90°. In effect, the direction of propagation is tangential to the 
surface of the plate, along the x-direction. The speed of sound is set at 343 m/s 
(1.3504 x 10 4 in. /s). The case of a grazing pressure distribution is no longer symmet- 
ric over the plate. The lower anti-symmetric modes, (1,2) and (2,1) must, therefore, 
be included in the modal transformation matrix. Similarly, the lowest membrane 
modes that were inactive under a normal load, may no longer be so. The plate model 
used in this case is of an 8 x 4 mesh with the first six bending modes [(1,1), (1,2), 
(2,1), (2,2), (1,3) and (3,1)] and four membrane modes. Thus, the total number of 
modes used for the grazing load model is R = 10. 

Response data are generated at 130db SPL. Apart from the center of the plate, 
deflection time-histories are recorded at every node along the x-axis centerline (a 
total of 15 non-zero nodal values). A linear simulation with grazing angle, (3 = 90°, 
with the six bending modes above is first carried out. The modal activity can now 
be compared against the modal results for (3 = 0°, normal load linear analysis. 
Given below are rms modal displacement, rms q r , values relative to rms q \ , the 
fundamental mode response, i.e., 


rms q r 
rms q\ 


r = 1,2,.. 


,6 
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(1,1) (1,2) (2,1) (2,2) (1,3) (3,1) 


(3 = 0° 1.0 0.0 0.0 0.0392 0.0948 0.0638 

/? = 90° 1.0 0.0261 0.1295 0.0350 0.0866 0.0279 


rmj q i(0°) 

rmj gi(90°) 


1.0099 


Table 9. Rms Modal ratios 


The contributions from each of these modes must be assessed with the individ- 
ual mode-shapes in mind (Figs.9 and 10). The participation from modes (1,2) and 
(2,1) when the load is modelled with angular incidence is expected. From the modal 
rms values above it is clear that response in mode (1,2) is very weak. Mode (2,1) is 
the strongest mode of response, apart from the fundamental mode. This mode, as 
seen in Fig. 9, has a trough and then a crest as one progresses along x-direction. Re- 
sponse over the plate in this shape, then, is a direct consequence of the left-to-right 
travelling wave assumption. The power spectrum for this grazing load response is 
shown in Fig.32. The spectrum at a node 3.75m. away from the leading edge is 
shown for both normal and grazing load models. This point was chosen as neither 
of the two additional modes, (1,2) and (2,1), will show up in the center deflection 
spectrum. Mode (1,2), is present in the response spectrum as a very small peak. A 
look at this mode in Fig. 9 will show that its amplitude along the x-axis centerline , 
is small. Compounded with the relatively low modal rms level for this mode, the 
modal contribution in the spectrum for that particular location is negligible. The 
same is true for mode (2,2). 

Deflection curves in Fig. 33 show the centerline rms deflections recorded at 
each node. They are plotted for both the normal and grazing load models. There 
is a slight change in the shape but overall response behavior shows no significant 
difference. Similarly, it was found that the PDF, power spectra and peak PDF 
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Fig. 32 Power Spectrum of Linear Deflection at x=3.75in. 
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distributions for maximum deflection showed little change. The power spectrum 
in the case of the grazing load model drops off a little more than for the normal 
load model, after about 500hz. The same trend was noticed in the strain response 
comparisons, too. None of these plot are displayed since the changes are rather 
insignificant. 

5.3.6 Variation of Other Parameters 

A few additional cases, involving parameters such as initial stresses from static 
compressive loads, iV 4 , quadrilateral elements, and plate thickness are studied and 
results are presented here. 

Plate Thickness 

The thickness parameter is studied by increasing the thickness of the baseline plate 
configuration ten-fold. The thickness of each layer is increased from 0.006 in. to 0.06 
in., resulting in total plate thickness, h = 0.48 in., with an a/h ratio of 31.25. The 
8x4 mesh with a total of six modes (4 bending and 2 membrane) is used. Linear 
and nonlinear simulations at 130db SPL are conducted. The centerline deflection 
values are cr//i — 0.002 for both the linear and nonlinear simulations. The strain 
data, too, were found to be the same for both the linear and nonlinear simulations. 
This implies that the responses of the plate are in the linear regime even at 130db 
SPL. The power spectral densities, as expected, were found to be identical in the 
two cases. As a matter of fact, although four bending modes were used, the modal 
response at the first mode seemed to be the strongest. 


Initial Stress 

The nonlinear random responses of the baseline panel under an initial stress dis- 
tribution are dealt with next. The initial stress distribution is accounted for using 
constant static compressive loads on all its edges. There is, thus, no initial moment 
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distribution in the plate. All other parameters of the simulation are the same as 
used in the previous example. Simulation is at 130db SPL using a total of six modes 
(4 bending and 2 membrane). The critical buckling static load for all four clamped 
edges under compression was computed at N x , N y = N cr = 20.23 lbs/in. There- 
fore) an edge load) = ]V C )- is used for this case. It may be recalled here that 
edge loads are set to be less than critical buckling load in the problem definition 
(see Eqs.2.3.4 and 2.3.9). 

The rms results for maximum deflection show the reduced stiffness in the plate 
due to this additional load. The value obtained is a/h = 1.88. The corresponding 
figure when edge loads are absent is 1.77. 

Quadrilateral Elements 

The Lagrangian element used in the finite element model in this study has so far 
been used only as a rectangular element. This element can actually be used, more 
generally, as a quadrilateral. This feature of the finite element model is demon- 
strated with an example. Fig. 34 displays the mesh structure used. The size of the 
mesh used was dictated by the necessity to generate quadrilateral shaped elements 
without any of them being too skewed. The mesh uses five elements along y-axis 
and 8,6 elements along the longer x-axis edges, respectively. The lower number of 
elements, along the y-axis makes the model stiffer along this direction, compared 
to the 8x8 mesh. The modal transformation matrix used, consisted of 4 bending 
modes and 2 membrane modes ( R = 6). In order to compare data from this model 
against the 8x8 baseline data, the four bending mode solution from the convergence 
studies (Fig. 12) is used here. Nonlinear simulation at 130db produced the following 
rms results 

a/h = 1.78 vs 1.78 for rectangular 8x8 mesh 

e rmJ = 497 fj.e vs 430 pe for rectangular 8x8 mesh 
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The deflection numbers above are compared at almost exactly the same loca- 
tion. The strain data are for elements of slightly different sizes. It is safe to say 
that most of the change in results is due to the change in mesh size. It was observed 
that the x-axis strain rms showed a similar degree of change. 
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Chapter 6 


CONCLUSIONS 

A Monte Carlo simulation technique to study the response of composite plates 
under combined acoustic and thermal loading has been implemented. Finite ele- 
ment equations in the physical degrees of freedom were reduced to a smaller set 
of modal equations of motion. An implicit numerical integration routine using a 
weighting function was encoded. A modified pressure distribution model was also 
generated. Response characteristics were obtained under realistic temperature dis- 
tribution assumptions. 

The modal decomposition of the large degree of freedom finite element model 
proved very effective in reducing the number of equations used for obtaining so- 
lutions. The formulation facilitated the inclusion of nonlinear stiffnesses without 
having to build them up at every time-step. Especially for a time-domain simu- 
lation procedure this technique is quite definitely a necessity. The computational 
time required otherwise will make the simulation method impractical. 

The numerical algorithm used appears to be a very useful tool [106]. The 
extension of the Wilson-0 method to the nonlinear dynamic equation has worked 
effectively. The algorithm is stable for fairly large nonlinearities at time-steps that 
are not too small, (At = 10~ 4 ). But the step-size cannot probably be reduced 
further since this is pretty close to the limit where amplitude decay and period 
elongation are small. The use of a nonlinear algebraic equation at each time-step 
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may be problem dependent. It was observed that for cases with moderate degrees 
of nonlinearity (low <r/h), solution of linearized equations at every time-step gave 
acceptable results. These results showed consistently higher rms deflection values. 
It is not likely, however, that strain results will be as accurate by using linearized 
equations within each time-step. In general, doing away with the iterations within 
time-steps, in effect becoming an explicit integration scheme, may be quite accept- 
able for systems with moderate nonlinearities. 

The modified load model was an effective tool in bringing out the complete 
spectrum of modal response. The implementation of this space-time correlated 
distribution did not add significantly to the computational cost. There was not 
much difference in the overall response characteristics, although the strain response 
shows a small degree of change. This is only true for the dimensions of the plate used 
in this study, relative to acoustic wavelengths. The grazing load model may show 
greater differences for a large or long plate. It will probably be more realistic to 
use a completely random (spatial and temporal) pressure distribution over the plate 
with this simulation scheme, since it can be implemented without much increase in 
computational effort. 

The results of simulations under combined acoustic and thermal loads have 
shown some response features not reported before. The influence of higher thermal 
buckling temperatures has been observed here, due to the unrestricted nature of 
the formulation. The data presented in Chapter 5, were at a high random load of 
130db. The same analysis at 90db, may have shown the behavior around buckling 
temperatures in a more dramatic fashion. The phenomenon of snap-through be- 
havior has also been simulated, although there is no way to approximately predict 
its onset. The surface variation used for the temperature increment gives a more 
realistic estimate of the temperature distribution in a plate. It would be interest- 
ing to record strain data at locations where the surface thermal gradients are the 
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steepest. Under a more severe thermal loading condition it is quite likely that the 
maximum strain location will no longer be close to the long-edge mid-point for this 
type of plate. 

The consideration of thermal moments induced by a temperature gradient 
through the plate thickness, does not seem to be much of a factor in the response 
behavior. Due to the direction of gradient used (normal to the heated surface is 
opposite to the dynamic load direction), there is probably a counteracting effect 
from the thermal moments. But the alteration in the response characteristics due 
to this gradient is still negligible. This is more so when the thermal load dominates 
the dynamic load, at low decibel levels [83]. 

A significant observation of this study has been the sensitivity of strain response 
characteristics to almost every parameter that influences the response behavior, 
including element type and mesh size used, membrane interpolation functions used, 
number of bending and membrane modes used, temperature distribution and load 
model used. In practical applications under such loading conditions as discussed in 
this work, it is the fatigue life performance of these structural components which 
is of critical importance. This needs accurate information on the strain response 
behavior of structures. Hence, it is important that strain response prediction models 
are accurate. 

In this context, the selection of mode-shapes to obtain converged solutions, is 
a critical part of the model definition process. In turn, it depends on the criterion 
chosen for modal convergence, namely, frequencies, deflection or strain. The num- 
ber of modes required, in turn, decides the mesh refinement procedure. Given the 
sensitivity of strain responses to all these factors, it is imperative that the conver- 
gence criterion, for both mesh refinement and modal convergence studies, be based 
on strain rather than deflections or natural frequencies. This tends to be the most 
time-consuming aspect of problem solution. The mesh required to provide accurate 
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strain results usually turns out to be large. In this study the eigenvector com- 
putation algorithms used were all full-system solution methods. The computation 
time for such solution processes for large degree-of-freedom systems, as we know, 
increases tremendously. It is therefore recommended for future applications that 
the eigenvalue-eigenvector computation algorithm used, be of the type that gener- 
ates a limited number of modes, say, one based on the Rayleigh quotient method. 
Once a good finite element modal transformation matrix has been obtained, the so- 
lution for any type of loading case for that particular laminate follows through with 
no restrictions. It is, therefore, recommended that future work focus on improved 
treatment of all these factors. 

Another aspect dealing with the particular type of laminate used in this study 
is with regard to boundary conditions. It has been observed that the presence of 
angle- ply layers removes the symmetry along the z, y-axes that a cross-ply laminate 
would have. This, therefore, precludes the use of quarter-plate models in the finite 
element mesh. The use of a quarter-plate model enforces symmetry conditions on 
all rotations and displacements, where no such symmetry exists. 

The simulation technique in combination with modal decomposition is proba- 
bly the best tool to build on, for analytical comparison with test data. It is also the 
best way to expand into coupled structural-acoustic prediction models. Its strength 
lies in its ability to model complicated problems with a minimum number of assump- 
tions and restrictions imposed. A vast library of time-histories can be used to give 
a wide range of statistical and spectral information that puts the response charac- 
teristics in clear perspective. In the process, one can observe behavioral phenomena 
not otherwise detected by other analytical approaches. Such a simulation, with a 
thorough modelling and statistical analysis of random responses under combined 
loading conditions is, possibly, the major contribution of this work. 
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Appendix A 


Finite Element Shape Functions 


The Lagrangian quadratic interpolation functions for a 9-node quadrilateral 
isoparametric element are listed below [135]. The variables £, rj axe natural coordi- 
nates whose range is — 1 < (£, 77 ) < 1 . They are defined such that 
9 9 


( x >y) = ( a ” > a i > a T \ a t* > a f v ) Nj 


1=1 


i=i 


are valid transformations. 


Ni 

N 2 
N 3 
N< 
N 5 
N 6 
N 7 
Ns 
N 9 


j(i - 0(i - ?)«i 

-jd -ew -v)i 
-j(i + 00 - n)(i 
-5(1 - 0(1 - -> 2 K 
(l - 0)0 - f) 
i(i +0(i -n 2 )i 
-j(l - 0(1 + '))('! 
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Appendix B 


Formation of Element Matrices 


Mass matrix 

This matrix is detailed first. The inertias, J, and shape functions, N t , have 
been defined earlier (Chapter 2 and Appendix A). 

elements 


[M] {a} = 
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h NiNj 
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d* 1 
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0 

0 

IzNiNj. 


CL j 


(i) 


where i, j = 1, . . . , 9; is the area over an element; [M c ] is a 45 x 45 element mass 
matrix consisting of sub-matrices [rn mm ], [m ww ] and [m**]. The element accelera- 
tion vector {a e } in turn consists of the corresponding sub- vectors a m , a w , anda^. 
The global matrix, [M] is of order Ndf- 
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hat 


Load Vector 

At a single time-step, say the n th time-step, the global load vector is given by 

f ° 1 

elements - I 0 I 

{/>-(()} = P „(() y ] N ' | dn ‘ 

elements 

= Pn (t) Y {*“} = P.(0{f} (2) 

where p n (<) is the magnitude of the pressure distribution at the n th time-step. The 
vector { F} is a global vector of length Ndf, representing unit pressure acting over 
the surface of the plate. The final form { P{t ) } is also of length Ndf- 

Note that of the 5 degrees of freedom per node, all terms for integration within 
an element (for the element load vector) except the one corresponding to the a w 
degree of freedom are set to zero. Load in the transverse direction is the only one 
being treated in this work. 


Linear Stiffness Matrices 

The terms in Eq.(2.4.10) generate the following submatrices at the element 
level. Indices p,r, s, t = 1, . . . , 9 and /, k — 1,...,18. 
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(7) 


are linear matrices corresponding to 5a™ rows. 
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where {/c om ] is the transpose of the matrix in Eq.(4), 
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and 
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( 10 ) 


all from w 0 , correspond to 8a™ rows. 



B.jt'pSKi d£T 


= f B„Ci,C jk dSl e a? 8af 
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( 11 ) 


where [fc^ m ] is the transpose of the matrix in Eq.(5), 



B, } e°8n t dn e 



B,jCuCjp dQ e a™ 8af 



(12) 


where [k^ 0 ] is the transpose of the matrix in Eq.(10), and 
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= *,?«!? *<■? 


(13) 


correspond to 


rows. Linear matrices from shear energy terms are 
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= / VA^B^BZ'dW a*S^ 
J n c 

= 

where [fcs u ' ,i ’] is the transpose of the matrix in Eq.(15), and 
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mnBZ r B” p dn' a p 
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(16) 


= ks™ w a™5a™ (17) 

Nonlinear Stiffness Matrices 

Other than one second-order term which is a cubic function of {fl ul }, all the 
other nonlinear energy expressions are first-order terms, quadratically dependent 
on {a™}- These are all listed below. 
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corresponds to 5a™ rows. The nonlinear strain tensor C- p3 has been made symmet- 
ric to be consistent with w, x w, y = w, y w, x . 



A tj e™5e f dCl e 
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(19) 


This is twice the transpose of kl 1 ^- . The second-order nonlinear function of {a“ 1 } 
is given by, 
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The other first-order functions of {a 11 '} are, 
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and 
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As in an earlier relation, Arirp, is half the transpose of fcl™° p . And 
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is half the transpose of This term corresponds to 8af rows. 

The constant terms in the strain energy are 
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Eqs.(3-29) can be put in a matrix form at the element level as follows 
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where all matrices and tensors are of order 45. 


+ < P ° 3 } ~ < P 


The linear stiffness matrix above is completely symmetric and so is the second- 
order nonlinear matrix. The submatrices corresponding to the {a m } and {a^} rows 
in the first-order matrix [kl] (12 and 32 positions) are transposes and half of the 
submatrices in the corresponding transpose positions of [fcl] (21 and 23 positions). 
The rest of the submatrices are symmetric. Thus, at the global assembled level we 
now have 


elements 
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E (l M 'H a '> + [C‘]{°'} + [A" + K 3t - K ATe ){a'} + [Al'JKH-'} 
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which results in the equation of motion 


[jlf]{a} + IC]{o} + [K + K‘ - K AT ][a} + [/n]{«}{«} 
+ [X2]{a}{a}{a} + {P'} - {P AT } 

= {P(<» 

as in Eq.(2.5.1). This system is of order N df . 
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Appendix C 


Numerical Integration of Modal Equations 


To proceed with deriving the nonlinear algebraic equations, the weighted in- 
tegral of the modal system of equations within one time-step is rewritten below. 
Note that p«(<) between the two time-steps, t n+J and has been replaced by p"(t) 
as in Eq.(3.1.8). The time subscripts n and (n + 1) are used as superscripts for 
subscripted variables. 


n + 1 f 

J W{t) j m, ] q j n+1 + c, } q^ +1 + ( k tJ + k 3 {j - k^ T )q^ +1 

t n 


+ 1 + Pi(t) + P? - Pi 1 j dr = 0 (1) 

Each one of the terms in the integrand will now be integrated in the manner of 
Eq.(3.1.5). The integration of the inertia term upon substitution of Eqs.(3.1.3) is 
as follows 
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( 2 ) 


= m,,(j; + aS J) "# 1 A() 


— m ij ( 9, + ^2 “7 / 


where d” and a^ n are related as in Eq.(3.1.9). The free parameter #i has been 
written in terms of the Wilson-# parameter, #w, as defined in Eqs.(3.1.7). 

In the same way, the damping term gives 
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Note that the usage of q } here includes the Wilson-# parameter with #i — #w- 
Similarly, the linear modal stiffness matrices together generate 
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(4) 


Here again g” includes = 6 w and 62 = 9\v in its expression. The {q} terms in 
Eqs.(2), (3) and (4) above, contribute to .AO, and the d n terms to A.,y. 

The quadratic and cubic tensor terms will be dealt with next. When using 
Eqs.(3.1.7) for the free parameters after integration, 9w and At occur together 
with like powers. This will be used in organizing the many terms as constant, 
linear, quadratic and cubic expressions of d n . 

The first-order nonlinear tensor term, thus results in 
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which are all constant terms contributing to ALO,, 


149 



r At 2 „ 9 A t\ 
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which is the quadratic term from Hot contributing to Alij k . 


In equation (5) 
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This property will be directly used in case of the second-order tensor whose 
integration gives 
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which are quadratic terms contributing to -41, jk, 

which is cubic in d n and gives A2,jjt/. 


The rest of the terms in (1) axe not dependent on r and hence remain un- 
changed. They add to .40, to form 

AO, = A0,((q},p,r°,p AT ) (8) 

The vector AOj remains unchanged with iterations and is only updated at every 
successive time-step. 

Therefore, grouping all the above terms as cubic, quadratic, linear and constant 
in d n , we arrive at 


A2ijki dfd^dj 1 + Al tjk d”d? + Aij d” + AO, = tf,(d) = 0 


(9) 
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